MSH3 ASO
Michael Flower
2024-10-31
Setup environment
# Set the working directory where all WGCNA functions are stored
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
getwd()## [1] "/Users/michaelflower/Library/Mobile Documents/com~apple~CloudDocs/Documents/ACL/Collaborations/Emma Bunting/2024.07.03 MSH3 ASO paper plots"Load packages.
# Required packages
packages <- c("plyr", "dplyr", "ggplot2", "data.table", "readxl", "tidyr", "tidyverse", "janitor", "patchwork",
"ggdark", "xlsx", "openxlsx", "ggpubr", "writexl", "stringr", "ggpmisc", "drc", "mgcv", "scales",
"rstatix", "kableExtra", "DT", "zoo", "emmeans", "lme4")
# # Install R packages if required
# install.packages(setdiff(packages, rownames(installed.packages())))
# # Install bioconductor packages if required
# BiocManager::install(setdiff(packages, rownames(installed.packages())))
# invisible(lapply(packages, function(x) library(x, character.only=TRUE)))
# Load
lapply(packages, library, character.only = TRUE)
rm(packages)
## Clear objects
# rm(list = ls()) # This would cause rmarkdown to error
# rm(list = ls()[!grepl("_env$", ls())]) # Remove all objects except environment importsSet variables
# Data path
data_path_spreadsheet <- "./data/MSH3aso_data.xlsx"
data_path_juliet <- "/Users/michaelflower/Library/Mobile Documents/com~apple~CloudDocs/Documents/ACL/Collaborations/Emma Bunting/2024.05.01 Emma MSH3 ASO FA reanalysis/results/02-juliet.RData"
# Output directory
out_dir <- "./results"
# Plot settings
element_text_size = 22
point_size = 3
hide_ns = TRUE
# Genes
gene_colour <- c("MSH2" = "red", # "#F8766D"
"MSH3" = "blue", # "#619CFF"
"MSH6" = "#00BA38")
gene_order <- c("MSH2", "MSH3", "MSH6")
# Treatments
treatment_rename <- c("Baseline" = "Baseline",
"Vehicle" = "Vehicle",
"MSH3aso-0.022uM" = "MSH3 ASO 0.022 µM",
"MSH3aso-0.26uM" = "MSH3 ASO 0.26 µM",
"MSH3aso-3uM" = "MSH3 ASO 3 µM",
"SCRaso-3uM" = "SCR ASO 3 µM",
"H2O2-500uM" = "H2O2 500 µM",
"SCRaso" = "SCR ASO",
"MSH3aso_1161149" = "MSH3 ASO 1161149",
"MSH3aso_1161173" = "MSH3 ASO 1161173",
"MSH3aso_1161329" = "MSH3 ASO 1161329")
treatment_order <- c("Baseline", "Vehicle", "MSH3aso-0.022uM", "MSH3aso-0.26uM",
"MSH3aso-3uM", "SCRaso-3uM", "SCRaso", "H2O2-500uM",
"MSH3aso_1161149", "MSH3aso_1161173", "MSH3aso_1161329")
treatment_colour <- c("Baseline" = "darkgrey",
"Vehicle" = "blue",
"MSH3aso-0.022uM" = "red",
"MSH3aso-0.26uM" = "purple",
"MSH3aso-3uM" = "orange",
"SCRaso-3uM" = "darkgreen",
"H2O2-500uM" = "cyan",
"SCRaso" = "darkgreen",
"MSH3aso_1161149" = "purple",
"MSH3aso_1161173" = "purple",
"MSH3aso_1161329" = "purple")
treatment_rename_subscript <- c("Baseline" = "Baseline",
"Vehicle" = "Vehicle",
"MSH3aso-0.022uM" = "MSH3~ASO~0.022~µM",
"MSH3aso-0.26uM" = "MSH3~ASO~0.26~µM",
"MSH3aso-3uM" = "MSH3~ASO~3~µM",
"SCRaso-3uM" = "SCR~ASO~3~µM",
"H2O2-500uM" = "H[2]*O[2]~500~µM")
# Genotypes
genotype_rename <- c("Unedited" = "Unedited",
"MSH3ko" = "MSH3-/-",
"FAN1ko" = "FAN1-/-")
genotype_order <- c("Unedited", "MSH3ko", "FAN1ko")
genotype_colour <- c("Unedited" = "blue",
"MSH3ko" = "darkred",
"FAN1ko" = "green")
# Ionis dose
treatment_ionis_rename <- c("Vehicle" = "Vehicle",
"3ug" = "MSH3 ASO 3 µg",
"10ug" = "MSH3 ASO 10 µg",
"30ug" = "MSH3 ASO 30 µg",
"SCRaso" = "SCR ASO")
treatment_ionis_order <- c("Vehicle", "3ug", "10ug", "30ug")
treatment_ionis_colour <- c("Vehicle" = "blue",
"3ug" = "red",
"10ug" = "purple",
"30ug" = "orange",
"SCRaso" = "darkgreen")
# Tissue
tissue_order <- c("Cortex", "Striatum", "Brainstem", "Spinal cord")
# Output formats
output_formats <- "svg" # c("png", "svg", "eps")
# Models to fit
my_models <- list(
"linear" = y ~ x,
"log" = y ~ log(x), # Can't use because instability rate goes negative at low MSH3 levels
"poly2" = y ~ poly(x, 2),
# "poly2raw" = y ~ poly(x, 2, raw = TRUE),
"poly3" = y ~ poly(x, 3),
# "poly3raw" = y ~ poly(x, 3, raw = TRUE),
"exponential" = I(log(y)) ~ x, # Exponential model # Can't use because instability rate goes negative at low MSH3 levels
# "logistic" = I(log(y / (1 - y))) ~ x, # Simple logistic model # Can't use because instability rate goes negative at low MSH3 levels
"exponential" = y ~ exp(x)
)
# Log increment
log_increment = 0.001
# p adjust method
padj_method = "BH" # "bonferroni"
# Errorbar metric
errorbar_metric = "SE" # CL, SDImport data
Spreadsheet
# Function to import all sheets
read_full_excel <- function(file_path = NULL){
# Imports the name of the sheets in the excel file and creates a vector
sheet_names <- readxl::excel_sheets(file_path)
# Create list and add each excel sheet as a dataframe by a loop
sheets <- list()
for(i in 1:length(sheet_names)){
sheets[[i]] <- readxl::read_xlsx(file_path, sheet = sheet_names[i])
}
# Add the sheet names to the list:
names(sheets) <- sheet_names
return(sheets)
}
# Import
MSH3aso_data <- read_full_excel(file_path = data_path_spreadsheet)Juliet
if (file.exists(file.path(out_dir, "juliet_minimal.RData"))) {
load(file.path(out_dir, "juliet_minimal.RData"))
} else {
# Import juliet data
if(!exists("juliet_env")) {
juliet_env <- new.env()
load(file = data_path_juliet,
envir = juliet_env)
ls(envir = juliet_env)
}
# Extract required objects
timecourse_stats <- juliet_env$timecourse_stats
timecourse_stats_manual <- juliet_env$timecourse_stats_manual
juliet_annotation <- juliet_env$annotation
# Extract slopes
my_stats <- timecourse_stats$predvar2$ii
my_experment_names <- setNames(names(my_stats), names(my_stats))
juliet_slopes <- lapply(my_experment_names, function(my_experiment_name,
my_stats) {
my_stats[[my_experiment_name]]$slopes_uncorrected %>%
dplyr::mutate(experiment_name = my_experiment_name) %>%
relocate(experiment_name)
},
my_stats = my_stats)
juliet_slopes <- rbindlist(juliet_slopes)
# Extract slopes from manually curated MSH3ko vs CRISPRwt experiment
my_slopes <- timecourse_stats_manual$predvar2$ii$slopes_uncorrected %>%
dplyr::mutate(experiment_name = "MSH3ko_CRISPRwt_manual") %>%
relocate(experiment_name)
juliet_slopes <- rbind(juliet_slopes, my_slopes)
# Save
save(juliet_slopes, juliet_annotation,
file = file.path(out_dir, "juliet_minimal.RData"))
}Model fitting function
Define a function to fit models.
# Vector of model names
model_names <- setNames(names(my_models), names(my_models))
# Function to plot fits
plot_model_function <- function(model_name,
models,
data,
xvar, xlabel,
yvar, ylabel,
group_var = NA, group_label = NA,
log_increment = 0.001,
plot_se = TRUE) {
# # Manual
# model_name = model_names[[2]]
# models = my_models
# data = plot_data
# xvar = "dose_uM"
# xlabel = "MSH3 ASO dose (µM)"
# yvar = "level"
# ylabel = "Relative protein level (%)"
# group_var = NA
# group_label = NA
# log_increment = 0.001
# plot_se = TRUE
# # Report
# print(model_name)
# Extract formula
my_model <- models[[model_name]]
# If formula contains log, use the pre-logged data with +0.001
my_model_str <- deparse(my_model)
contains_log <- grepl("log", my_model_str)
if (contains_log) {
data <- data %>%
dplyr::mutate(!!sym(xvar) := !!sym(xvar) + log_increment)
}
# Plot
if (is.na(group_var)) {
my_plot <-
ggplot(data, aes(x = !!sym(xvar), y = !!sym(yvar))) +
# geom_point(size = point_size) +
geom_point() +
stat_smooth(method = lm,
formula = my_model,
fullrange = TRUE, se = plot_se) +
stat_poly_eq(formula = my_model,
aes(label = paste(after_stat(eq.label),
after_stat(rr.label),
after_stat(p.value.label),
sep = "~~~")),
parse = TRUE, coef.digits = 3, f.digits = 3, p.digits = 3,
rr.digits = 3, size = 3) +
labs(title = model_name,
subtitle = c(my_model),
x = xlabel,
y = ylabel) +
theme_minimal()
} else {
my_plot <-
ggplot(data, aes(x = !!sym(xvar), y = !!sym(yvar), colour = !!sym(group_var), fill = !!sym(group_var))) +
geom_point() +
stat_smooth(method = lm,
formula = my_model,
fullrange = TRUE, se = plot_se) +
stat_poly_eq(formula = my_model,
aes(label = paste(after_stat(eq.label),
after_stat(rr.label),
after_stat(p.value.label),
sep = "~~~")),
parse = TRUE, coef.digits = 3, f.digits = 3, p.digits = 3,
rr.digits = 3, size = 3) +
labs(title = model_name,
subtitle = c(my_model),
x = xlabel,
y = ylabel,
colour = group_label,
fill = group_label) +
theme_minimal()
}
# Convert xy notation to real variables for lm function
left <- gsub("\\by\\b", yvar, as.character(my_model[2]))
right <- gsub("\\bx\\b", xvar, as.character(my_model[3])) # This pattern \\bx\\b ensures that "x" is replaced only when it stands alone as a word. \\b denotes a word boundary.
# Linear model
if (is.na(group_var)) {
my_model_long <- sprintf("%s ~ %s", left, right) # https://stackoverflow.com/questions/26381410/edit-and-reuse-the-formula-part-of-the-call-for-a-model-in-r
my_model_long <- as.formula(my_model_long)
my_lm <- lm(my_model_long, data = data)
r2 <- summary(my_lm)$r.squared
} else {
my_model_long <- sprintf("%s ~ %s * %s", left, right, group_var) # https://stackoverflow.com/questions/26381410/edit-and-reuse-the-formula-part-of-the-call-for-a-model-in-r
my_model_long <- as.formula(my_model_long)
my_lm <- lm(my_model_long, data = data)
r2 <- summary(my_lm)$r.squared
}
# Output
out <- mget(c("my_plot", "my_model", "r2"))
return(out)
}Titration western MSH3 1wk
Fit models
## [1] "titration_western_MSH3_1wk"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!dose_uM,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("gene", "rep")) %>%
dplyr::mutate(dose_uM = as.numeric(dose_uM),
level = as.numeric(level),
dose_uM_forlog = dose_uM + log_increment) %>% # Add increment (+0.001) to log the data
dplyr::filter(!is.na(level))
# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
models = my_models,
data = plot_data,
xvar = "dose_uM",
xlabel = "MSH3 ASO dose (µM)",
yvar = "level",
ylabel = "Relative protein level (%)",
group_var = NA,
group_label = NA,
log_increment = log_increment,
plot_se = TRUE)
# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
my_fit$my_plot
})
wrap_plots(my_plots) +
plot_annotation(title = "Regression models",
theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))
# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
r2 = unlist(my_r2)) %>%
tibble::rownames_to_column("model_number") %>%
arrange(-r2)
# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
kable_styling(bootstrap_options = "striped")| model_number | model_name | r2 |
|---|---|---|
| 2 | log | 0.8776671 |
| 4 | poly3 | 0.8705157 |
| 5 | exponential | 0.7160194 |
| 6 | exponential | 0.7160194 |
| 3 | poly2 | 0.6906890 |
| 1 | linear | 0.5987664 |
Functions
# Helper function to calculate the effective dose for a given response level using interpolation
ed_interpolation_helper <- function(data, response_var, pred_var, formula, model_function, response_level) {
if (model_function == "lm") {
my_model <- lm(formula, data = data)
} else if (model_function == "drm") {
my_model <- drm(formula, data = data, fct = LL.4())
}
predicted_values <- setNames(data.frame(seq(0, max(data[[pred_var]]), length.out = 1000)), pred_var)
predictions <- predict(my_model, newdata = predicted_values, interval = "confidence")
predicted_values <- cbind(predicted_values, predictions)
colnames(predicted_values) <- c(pred_var, response_var, "lower_ci", "upper_ci")
ed <- approx(predicted_values[, response_var], predicted_values[[pred_var]], xout = response_level)$y
return(list(predicted_values = predicted_values,
ed = ed))
}
# Effective dose bootstrap function to calculate confidence intervals
ed_bootstrap_helper <- function(data, response_var, pred_var, formula, model_function, response_level, n_bootstrap = 1000, confidence_level = 0.95) {
# Resample data and compute ED for each resample
ed_bs <- replicate(n_bootstrap, {
data_resample <- data[sample(nrow(data), replace = TRUE), ]
my_ed <- ed_interpolation_helper(data = data_resample,
response_var = response_var,
pred_var = pred_var,
formula = formula,
model_function = model_function,
response_level = response_level)
return(my_ed$ed)
})
# Calculate confidence intervals
alpha <- 1 - confidence_level
ci_lower = quantile(na.omit(ed_bs), probs = alpha / 2)
ci_upper = quantile(na.omit(ed_bs), probs = 1 - (alpha / 2))
# Output
return(list(ci_lower = as.numeric(ci_lower),
ci_upper = as.numeric(ci_upper)))
}
# Bring the effective dose (ED) interpolation and confidence interval functions together
ed_interpolation <- function(data, response_var, pred_var, formula, model_function, response_level, n_bootstrap = 1000, confidence_level = 0.95) {
# Calculate effective dose
ed_result <- ed_interpolation_helper(data = data,
response_var = response_var,
pred_var = pred_var,
formula = formula,
model_function = model_function,
response_level = response_level)
# Confidence intervals
ed_ci <- ed_bootstrap_helper(data = data,
response_var = response_var,
pred_var = pred_var,
formula = formula,
model_function = model_function,
response_level = response_level,
n_bootstrap = n_bootstrap,
confidence_level = confidence_level)
# Output
return(list(ed = ed_result$ed,
ci_lower = ed_ci$ci_lower,
ci_upper = ed_ci$ci_upper,
predicted_values = ed_result$predicted_values))
}IC50
# Dose response model
my_drm <- drm(level ~ dose_uM, data = plot_data, fct = LL.4())
# Extract IC50 value
ic50 <- ED(my_drm, 50)##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 1.7779 10.8055ic50_estimate <- ic50[[1]]
ic50_se <- ic50[[2]]
ic50_upper <- ic50_estimate + 1.96 * ic50_se
ic50_lower <- ic50_estimate - 1.96 * ic50_se
# Plot drm model
plot(my_drm, main = "Dose-response curve using drc")
abline(h = 50, col = "blue", lty = 2)
abline(v = ic50_estimate, col = "blue", lty = 2)
# Calculate IC50 by interpolation
ic50_interpolation <- ed_interpolation(data = plot_data,
response_var = "level",
pred_var = "dose_uM_forlog",
formula = as.formula("level ~ log(dose_uM_forlog)"),
model_function = "lm",
response_level = 50,
n_bootstrap = 1000,
confidence_level = 0.95)
# Readjust for the log increment
ic50_interpolation$ed <- ic50_interpolation$ed - log_increment
ic50_interpolation$ci_lower <- ic50_interpolation$ci_lower - log_increment
ic50_interpolation$ci_upper <- ic50_interpolation$ci_upper - log_increment
# Plot linear model
plot(level ~ dose_uM_forlog,
data = ic50_interpolation$predicted_values,
main = "Dose-response curve using interpolation")
abline(h = 50, col = "blue", lty = 2)
abline(v = ic50_interpolation$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation$ci_upper, col = "lightblue", lty = 2)
# Table of IC50
ic50_table <- data.frame(method = c("drc", "lm_interpolation"),
IC50 = c(ic50_estimate, ic50_interpolation$ed),
ci_lower = c(ic50_lower, ic50_interpolation$ci_lower),
ci_upper = c(ic50_upper, ic50_interpolation$ci_upper))
# Display IC50 table
kbl(ic50_table, caption = "IC50") %>%
kable_styling(bootstrap_options = "striped")| method | IC50 | ci_lower | ci_upper |
|---|---|---|---|
| drc | 1.7779260 | -19.4008695 | 22.9567214 |
| lm_interpolation | 0.2108101 | 0.1296216 | 0.3853176 |
Linear dose axis
# Plot
my_plot <-
ggplot(plot_data, aes(x = dose_uM_forlog, y = level)) +
geom_point(aes(colour = gene), size = point_size) +
stat_smooth(aes(colour = gene, fill = gene),
method = lm, formula = y ~ log(x), fullrange = TRUE,
size = 1, alpha = 0.1) +
geom_segment(aes(x = -Inf, xend = ic50_interpolation$ed,
y = 50, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed,
y = -Inf, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
ymin = -Inf, ymax = 50),
fill = "darkgrey", alpha = 0.03, colour = NA) +
scale_colour_manual(values = gene_colour) +
scale_fill_manual(values = gene_colour) +
scale_x_continuous(breaks = pretty_breaks()) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 ASO dose (µM)",
y = "Relative MSH3 protein level (%)",
colour = "Protein",
fill = "Protein") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx.", output_format)),
plot = my_plot)
}## Saving 7 x 5 in image
Log dose axis
# Plot
my_plot <-
ggplot(plot_data, aes(x = dose_uM_forlog, y = level, colour = gene, fill = gene)) +
geom_point(size = point_size) +
stat_smooth(method = lm,
formula = y ~ x, fullrange = TRUE,
size = 1, alpha = 0.1) +
geom_segment(aes(x = 0.001, xend = ic50_interpolation$ed,
y = 50, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed,
y = -Inf, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
ymin = -Inf, ymax = 50),
fill = "darkgrey", alpha = 0.03, colour = NA) +
scale_colour_manual(values = gene_colour) +
scale_fill_manual(values = gene_colour) +
scale_x_continuous(trans = "log2",
breaks = c(0.001, 0.01, 0.1, 1, 2, 3),
labels = function(x) {
ifelse(x < 1, round(x, digits = 3), round(x, digits = 0))
}) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 ASO dose (µM)",
y = "Relative MSH3 protein level (%)",
colour = "Protein",
fill = "Protein") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx.", output_format)),
plot = my_plot)
}## Saving 7 x 5 in image
Titration western MSH2 and MSH6 1wk
Fit models
## [1] "titration_western_MSH2_MSH6_1wk"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!dose_uM,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("gene", "rep")) %>%
dplyr::mutate(dose_uM = as.numeric(dose_uM),
level = as.numeric(level),
dose_uM_forlog = dose_uM + log_increment) %>% # Add increment (+0.001) to log the data
dplyr::filter(!is.na(level))
# Relevel variables
plot_data <- plot_data %>%
dplyr::mutate(gene = fct_relevel(gene, intersect(gene_order, unique(plot_data$gene))))
# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
models = my_models,
data = plot_data,
xvar = "dose_uM",
xlabel = "MSH3 ASO dose (µM)",
yvar = "level",
ylabel = "Relative protein level (%)",
group_var = "gene",
group_label = "Gene",
log_increment = log_increment,
plot_se = TRUE)
# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
my_fit$my_plot
})
wrap_plots(my_plots) +
plot_annotation(title = "Regression models",
theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))
# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
r2 = unlist(my_r2)) %>%
tibble::rownames_to_column("model_number") %>%
arrange(-r2)
# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
kable_styling(bootstrap_options = "striped")| model_number | model_name | r2 |
|---|---|---|
| 4 | poly3 | 0.1953677 |
| 3 | poly2 | 0.1478138 |
| 2 | log | 0.1371533 |
| 5 | exponential | 0.1025985 |
| 6 | exponential | 0.1025985 |
| 1 | linear | 0.1009314 |
Linear dose axis
# Plot
my_plot <-
ggplot(plot_data, aes(x = dose_uM, y = level, colour = gene, fill = gene)) +
geom_point(size = point_size) +
stat_smooth(method = lm,
formula = y ~ x,
fullrange = TRUE, alpha = 0.2) +
scale_colour_manual(values = gene_colour) +
scale_fill_manual(values = gene_colour) +
scale_x_continuous(breaks = pretty_breaks()) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 ASO dose (µM)",
y = "Relative protein level (%)",
colour = "Protein",
fill = "Protein") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx.", output_format)),
plot = my_plot)
}## Saving 7 x 5 in image
Log dose axis
# Plot
my_plot <-
ggplot(plot_data, aes(x = dose_uM, y = level, colour = gene, fill = gene)) +
geom_point(size = point_size) +
stat_smooth(method = lm,
formula = y ~ x,
fullrange = TRUE, alpha = 0.2) +
scale_colour_manual(values = gene_colour) +
scale_fill_manual(values = gene_colour) +
scale_x_continuous(trans = "log2",
breaks = c(0.001, 0.01, 0.1, 1, 2, 3),
labels = function(x) {
ifelse(x < 1, round(x, digits = 3), round(x, digits = 0))
}) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 ASO dose (µM)",
y = "Relative protein level (%)",
colour = "Protein",
fill = "Protein") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx.", output_format)),
plot = my_plot)
}## Saving 7 x 5 in image
Titration western MSH3 9wk
Fit models
## [1] "titration_western_MSH3_9wk"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!dose_uM,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("gene", "rep")) %>%
dplyr::mutate(dose_uM = as.numeric(dose_uM),
level = as.numeric(level),
dose_uM_forlog = dose_uM + log_increment) %>% # Add increment (+0.001) to log the data
dplyr::filter(!is.na(level))
# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
models = my_models,
data = plot_data,
xvar = "dose_uM",
xlabel = "MSH3 ASO dose (µM)",
yvar = "level",
ylabel = "Relative protein level (%)",
group_var = NA,
group_label = NA,
log_increment = log_increment,
plot_se = TRUE)
# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
my_fit$my_plot
})
wrap_plots(my_plots) +
plot_annotation(title = "Regression models",
theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))
# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
r2 = unlist(my_r2)) %>%
tibble::rownames_to_column("model_number") %>%
arrange(-r2)
# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
kable_styling(bootstrap_options = "striped")| model_number | model_name | r2 |
|---|---|---|
| 4 | poly3 | 0.9217185 |
| 2 | log | 0.9204592 |
| 3 | poly2 | 0.8405180 |
| 5 | exponential | 0.7286988 |
| 6 | exponential | 0.7286988 |
| 1 | linear | 0.5735873 |
IC50
# Dose response model
my_drm <- drm(level ~ dose_uM, data = plot_data, fct = LL.4())
# Extract IC50 value
ic50 <- ED(my_drm, 50)##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 0.13746 0.35855ic50_estimate <- ic50[[1]]
ic50_se <- ic50[[2]]
ic50_upper <- ic50_estimate + 1.96 * ic50_se
ic50_lower <- ic50_estimate - 1.96 * ic50_se
# Plot drm model
plot(my_drm, main = "Dose-response curve using drc")
abline(h = 50, col = "green", lty = 2)
abline(v = ic50_estimate, col = "blue", lty = 2)
# Calculate IC50 by interpolation
ic50_interpolation <- ed_interpolation(data = plot_data,
response_var = "level",
pred_var = "dose_uM_forlog",
formula = as.formula("level ~ log(dose_uM_forlog)"),
model_function = "lm",
response_level = 50,
n_bootstrap = 1000,
confidence_level = 0.95)
# Readjust for the log increment
ic50_interpolation$ed <- ic50_interpolation$ed - log_increment
ic50_interpolation$ci_lower <- ic50_interpolation$ci_lower - log_increment
ic50_interpolation$ci_upper <- ic50_interpolation$ci_upper - log_increment
# Plot linear model
plot(level ~ dose_uM_forlog,
data = ic50_interpolation$predicted_values,
main = "Dose-response curve using interpolation")
abline(h = 50, col = "blue", lty = 2)
abline(v = ic50_interpolation$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation$ci_upper, col = "lightblue", lty = 2)
# Table of IC50
ic50_table <- data.frame(method = c("drc", "lm_interpolation"),
IC50 = c(ic50_estimate, ic50_interpolation$ed),
ci_lower = c(ic50_lower, ic50_interpolation$ci_lower),
ci_upper = c(ic50_upper, ic50_interpolation$ci_upper))
# Display IC50 table
kbl(ic50_table, caption = "IC50") %>%
kable_styling(bootstrap_options = "striped")| method | IC50 | ci_lower | ci_upper |
|---|---|---|---|
| drc | 0.1374605 | -0.5652906 | 0.8402115 |
| lm_interpolation | 0.0915830 | 0.0436033 | 0.1893050 |
Linear dose axis
# Plot
my_plot <-
ggplot(plot_data, aes(x = dose_uM_forlog, y = level, colour = gene, fill = gene)) +
geom_point(size = point_size) +
stat_smooth(method = lm,
formula = y ~ log(x), fullrange = TRUE,
size = 1, alpha = 0.1) +
geom_segment(aes(x = -Inf, xend = ic50_interpolation$ed,
y = 50, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed,
y = -Inf, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
ymin = -Inf, ymax = 50),
fill = "darkgrey", alpha = 0.03, colour = NA) +
# geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
# ymin = 0, ymax = Inf), fill = "grey", alpha = 0.05, colour = NA) +
# geom_hline(yintercept = 50, colour = "black", linetype = "dashed") +
# geom_vline(xintercept = ic50_interpolation$ed, colour = "black", linetype = "dashed") +
scale_colour_manual(values = gene_colour) +
scale_fill_manual(values = gene_colour) +
scale_x_continuous(breaks = pretty_breaks()) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 ASO dose (µM)",
y = "Relative MSH3 protein level (%)",
colour = "Protein",
fill = "Protein") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_linearx.", output_format)),
plot = my_plot)
}## Saving 7 x 5 in image
Log dose axis
# Plot
my_plot <-
ggplot(plot_data, aes(x = dose_uM_forlog, y = level, colour = gene, fill = gene)) +
geom_point(size = point_size) +
stat_smooth(method = lm,
formula = y ~ x, fullrange = TRUE,
size = 1, alpha = 0.1) +
geom_segment(aes(x = 0.001, xend = ic50_interpolation$ed,
y = 50, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_segment(aes(x = ic50_interpolation$ed, xend = ic50_interpolation$ed,
y = -Inf, yend = 50),
colour = "black", linetype = "dashed", size = 0.5) +
geom_rect(aes(xmin = ic50_interpolation$ci_lower, xmax = ic50_interpolation$ci_upper,
ymin = -Inf, ymax = 50),
fill = "darkgrey", alpha = 0.03, colour = NA) +
scale_colour_manual(values = gene_colour) +
scale_fill_manual(values = gene_colour) +
scale_x_continuous(trans = "log2",
breaks = c(0.001, 0.01, 0.1, 1, 2, 3),
labels = function(x) {
ifelse(x < 1, round(x, digits = 3), round(x, digits = 0))
}) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 ASO dose (µM)",
y = "Relative MSH3 protein level (%)",
colour = "Protein",
fill = "Protein") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# I'M HERE!
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_logx.", output_format)),
plot = my_plot)
}## Saving 7 x 5 in image
Duration response
Fit models
## [1] "duration_response"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!treatment_duration_weeks,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
dplyr::mutate(treatment_duration_weeks = as.numeric(treatment_duration_weeks),
level = as.numeric(level) * 100) %>%
dplyr::filter(!is.na(level))
# Relevel variables
plot_data <- plot_data %>%
dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))))
# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
models = my_models,
data = plot_data,
xvar = "treatment_duration_weeks",
xlabel = "Week",
yvar = "level",
ylabel = "Relative expression (%)",
group_var = "treatment",
group_label = "Treatment",
log_increment = log_increment,
plot_se = TRUE)
# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
my_fit$my_plot
})
wrap_plots(my_plots) +
plot_annotation(title = "Regression models",
theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))
# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
r2 = unlist(my_r2)) %>%
tibble::rownames_to_column("model_number") %>%
arrange(-r2)
# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
kable_styling(bootstrap_options = "striped")| model_number | model_name | r2 |
|---|---|---|
| 4 | poly3 | 0.7813014 |
| 3 | poly2 | 0.7771263 |
| 2 | log | 0.7704015 |
| 1 | linear | 0.7668544 |
| 5 | exponential | 0.7420702 |
| 6 | exponential | 0.7420702 |
Plot
# Plot
my_plot <-
ggplot(plot_data, aes(x = treatment_duration_weeks, y = level, colour = treatment, fill = treatment)) +
geom_point(size = point_size) +
stat_smooth(method = lm,
formula = y ~ log(x),
fullrange = FALSE, alpha = 0.2) +
scale_colour_manual(values = treatment_colour,
labels = treatment_rename) +
scale_fill_manual(values = treatment_colour,
labels = treatment_rename) +
scale_x_continuous(breaks = c(0, 3, 6, 9, 12, 15)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Weeks",
y = "Relative expression (%)",
colour = "Treatment",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
plot = my_plot)
}## Saving 10 x 7 in image
Statistics
Descriptive statistics
Mean, SD and SEM MSH3 level at each timepoint, for each treatment.
# Descriptive statistics
descriptive_stats <- plot_data %>%
group_by(treatment_duration_weeks, treatment) %>%
summarise(mean = mean(level),
sd = sd(level),
sem = sd(level)/sqrt(n())) %>%
ungroup()## `summarise()` has grouped output by 'treatment_duration_weeks'. You can override using the `.groups` argument.# Display table
kbl(descriptive_stats, caption = "Descriptive statistics") %>%
kable_styling(bootstrap_options = "striped")| treatment_duration_weeks | treatment | mean | sd | sem |
|---|---|---|---|---|
| 3 | Vehicle | 100.000000 | 0.000000 | 0.0000000 |
| 3 | MSH3aso-0.022uM | 103.753600 | 42.175955 | 24.3502992 |
| 3 | MSH3aso-0.26uM | 19.315167 | 3.901393 | 2.2524702 |
| 3 | MSH3aso-3uM | 5.976033 | 7.204422 | 4.1594753 |
| 3 | SCRaso-3uM | 190.146533 | 99.644166 | 57.5295859 |
| 9 | Vehicle | 100.000000 | 0.000000 | 0.0000000 |
| 9 | MSH3aso-0.022uM | 50.742850 | 39.062771 | 27.6215500 |
| 9 | MSH3aso-0.26uM | 47.476467 | 35.675978 | 20.5975358 |
| 9 | MSH3aso-3uM | 1.235200 | NA | NA |
| 9 | SCRaso-3uM | 159.393400 | 51.249458 | 29.5888885 |
| 12 | Vehicle | 100.000000 | 0.000000 | 0.0000000 |
| 12 | MSH3aso-0.022uM | 34.276467 | 9.674190 | 5.5853963 |
| 12 | MSH3aso-0.26uM | 20.043367 | 23.866658 | 13.7794212 |
| 12 | MSH3aso-3uM | 0.678400 | 1.112659 | 0.6423939 |
| 12 | SCRaso-3uM | 150.230967 | 36.021422 | 20.7969779 |
| 15 | Vehicle | 100.000000 | 0.000000 | 0.0000000 |
| 15 | MSH3aso-0.022uM | 46.574467 | 55.495055 | 32.0400852 |
| 15 | MSH3aso-0.26uM | 23.504833 | 36.844000 | 21.2718932 |
| 15 | MSH3aso-3uM | 3.307650 | 4.641237 | 3.2818500 |
| 15 | SCRaso-3uM | 132.536300 | 54.266488 | 31.3307716 |
Group mean
This fits a repeated measures linear mixed-effects model ANOVA,
lmer(level ~ treatment * treatment_duration_weeks + (1|rep)),
to control for changes over time and different biological replicates. It
calculates the mean MSH3 levels for each treatment using estimated
marginal means and performs posthoc pairwise comparisons. The results
are displayed in a table, and the group means are visualised in a
plot.
# Repeated measures linear mixed-effects model ANOVA
my_lm <- lmer(level ~ treatment * treatment_duration_weeks + (1|rep), data = plot_data)
anova(my_lm)# Number of timecourses contributing to each group
sample_count <- plot_data %>%
distinct(treatment_duration_weeks, treatment) %>%
group_by(treatment) %>%
dplyr::summarise(n_timepoints = n()) %>%
ungroup()
# Posthoc test calculating and comparing group averages
my_emm <- emmeans(my_lm, pairwise ~ treatment | treatment_duration_weeks) # Post-hoc tests of group averages
# Extract group means
group_mean <- as.data.frame(summary(my_emm)$emmeans) %>%
left_join(sample_count, by = join_by(treatment)) %>%
dplyr::mutate(SD = SE * sqrt(n_timepoints)) %>%
relocate(SD, .before = "SE")
# Display group means
kbl(group_mean, caption = "Mean MSH3 level for each treatment group") %>%
kable_styling(bootstrap_options = "striped")| treatment | treatment_duration_weeks | emmean | SD | SE | df | lower.CL | upper.CL | n_timepoints |
|---|---|---|---|---|---|---|---|---|
| Vehicle | 9.7 | 100.000000 | 18.08702 | 9.043511 | 14.36387 | 80.649594 | 119.35041 | 4 |
| MSH3aso-0.022uM | 9.7 | 61.291940 | 21.64865 | 10.824327 | 21.37144 | 38.805318 | 83.77856 | 4 |
| MSH3aso-0.26uM | 9.7 | 28.036304 | 20.63006 | 10.315032 | 19.98939 | 6.518792 | 49.55382 | 4 |
| MSH3aso-3uM | 9.7 | 3.497054 | 23.82746 | 11.913730 | 26.05924 | -20.989259 | 27.98337 | 4 |
| SCRaso-3uM | 9.7 | 158.764153 | 20.63006 | 10.315032 | 19.98939 | 137.246642 | 180.28166 | 4 |
# Pairwise comparison of treatment groups
group_mean_compare <- as.data.frame(summary(my_emm)$contrasts) %>%
dplyr::rename(p.adj = p.value) %>% # emmeans already adjusts for multiple comparisons using tukey
tidyr::separate(contrast, c("group1", "group2"), sep=" - ") %>%
dplyr::mutate(p.signif = symnum(p.adj,
corr = FALSE,
na = FALSE,
legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", "")),
group1 = gsub("^\\(|\\)$", "", group1), # Remove brackets
group2 = gsub("^\\(|\\)$", "", group2)) # Remove brackets
# Hide non-significant results
if (hide_ns) {
group_mean_compare_plot <- group_mean_compare %>%
filter(as.numeric(p.adj) < 0.05)
} else {
group_mean_compare_plot <- pwc_result
}
# Display pairwise comparisons
kbl(group_mean_compare, caption = "Pairwise comparison of treatment group means") %>%
kable_styling(bootstrap_options = "striped")| group1 | group2 | treatment_duration_weeks | estimate | SE | df | t.ratio | p.adj | p.signif |
|---|---|---|---|---|---|---|---|---|
| Vehicle | MSH3aso-0.022uM | 9.7 | 38.70806 | 12.91980 | 48.38705 | 2.996025 | 0.0333704 |
|
| Vehicle | MSH3aso-0.26uM | 9.7 | 71.96370 | 12.49621 | 47.78943 | 5.758844 | 0.0000057 | *** |
| Vehicle | MSH3aso-3uM | 9.7 | 96.50295 | 13.84530 | 48.69702 | 6.970089 | 0.0000001 | *** |
| Vehicle | SCRaso-3uM | 9.7 | -58.76415 | 12.49621 | 47.78943 | -4.702560 | 0.0002079 | *** |
| MSH3aso-0.022uM | MSH3aso-0.26uM | 9.7 | 33.25564 | 13.57150 | 47.29441 | 2.450402 | 0.1198478 | |
| MSH3aso-0.022uM | MSH3aso-3uM | 9.7 | 57.79489 | 14.64059 | 47.28273 | 3.947579 | 0.0023225 | ** |
| MSH3aso-0.022uM | SCRaso-3uM | 9.7 | -97.47221 | 13.57150 | 47.29441 | -7.182124 | 0.0000000 | *** |
| MSH3aso-0.26uM | MSH3aso-3uM | 9.7 | 24.53925 | 14.39221 | 47.44504 | 1.705037 | 0.4407232 | |
| MSH3aso-0.26uM | SCRaso-3uM | 9.7 | -130.72785 | 13.24338 | 47.09175 | -9.871183 | 0.0000000 | *** |
| MSH3aso-3uM | SCRaso-3uM | 9.7 | -155.26710 | 14.39221 | 47.44504 | -10.788273 | 0.0000000 | *** |
# Plot
my_plot <-
ggplot(group_mean, aes(x = treatment, y = emmean)) +
geom_bar(stat = "identity",
colour = "black",
# aes(fill = treatment), alpha = 0.8) +
aes(fill = treatment)) +
geom_jitter(data = plot_data,
aes(x = treatment, y = level, fill = treatment),
height = 0, width = 0.1,
size = 3, shape = 21, colour = "black", alpha = 0.8) +
geom_errorbar(aes(ymin = case_when(errorbar_metric == "CL" ~ lower.CL,
errorbar_metric == "SE" ~ emmean - SE,
errorbar_metric == "SD" ~ emmean - SD,
TRUE ~ NA),
ymax = case_when(errorbar_metric == "CL" ~ upper.CL,
errorbar_metric == "SE" ~ emmean + SE,
errorbar_metric == "SD" ~ emmean + SD,
TRUE ~ NA)),
width = 0.2) +
stat_pvalue_manual(group_mean_compare_plot,
label = "{p.signif}",
y.position = ifelse(errorbar_metric == "CL",
1.1 * max(group_mean$upper.CL),
1.1 * (max(group_mean$emmean) + max(group_mean$SE))),
step.increase = 0.1,
vjust = 0.75,
hide.ns = hide_ns,
size = 6) +
scale_x_discrete(labels = treatment_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
scale_colour_manual(values = treatment_colour) +
scale_fill_manual(values = treatment_colour,) +
labs(x = "Treatment",
y = "Relative MSH3 expression (%)",
colour = "Treatment",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
# Dispay plot
print(my_plot)
Group slope
This calculates and compares the slopes of change in MSH3 levels for each treatment. The slopes are estimated using the emtrends function from the emmeans package, which fits a repeated measures linear mixed-effects model ANOVA. Posthoc pairwise comparisons are performed to compare the slopes between treatment groups. The results are displayed in a table, and the slopes are visualised in a plot.
# Posthoc test calculating and comparing group slopes
my_emt <- emtrends(my_lm, ~ treatment, var = "treatment_duration_weeks")
# Extract group slopes
group_slope <- as.data.frame(my_emt)
# Display group slopes
kbl(group_slope, caption = "Slope of change in MSH3 level for each treatment group") %>%
kable_styling(bootstrap_options = "striped")| treatment | treatment_duration_weeks.trend | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| Vehicle | 0.0000000 | 1.827644 | 47.09175 | -3.676557 | 3.6765570 |
| MSH3aso-0.022uM | -5.0192411 | 2.158652 | 48.92396 | -9.357385 | -0.6810973 |
| MSH3aso-0.26uM | 0.3708782 | 2.160993 | 49.08288 | -3.971614 | 4.7133708 |
| MSH3aso-3uM | -0.0562575 | 2.289665 | 48.28564 | -4.659234 | 4.5467194 |
| SCRaso-3uM | -4.3492758 | 2.160993 | 49.08288 | -8.691768 | -0.0067832 |
# Pairwise comparisons
my_pwc <- pairs(my_emt)
group_slope_compare <- data.frame(my_pwc) %>%
dplyr::rename(p.adj = p.value) %>% # emtrends already adjusts for multiple comparisons using tukey
tidyr::separate(contrast, c("group1", "group2"), sep=" - ") %>%
dplyr::mutate(p.signif = symnum(p.adj,
corr = FALSE,
na = FALSE,
legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", "")))
# Display pairwise comparisons
kbl(group_slope_compare, caption = "Pairwise comparison of treatment group slopes") %>%
kable_styling(bootstrap_options = "striped")| group1 | group2 | estimate | SE | df | t.ratio | p.adj | p.signif |
|---|---|---|---|---|---|---|---|
| Vehicle | (MSH3aso-0.022uM) | 5.0192411 | 2.828438 | 48.24414 | 1.7745629 | 0.3998857 | |
| Vehicle | (MSH3aso-0.26uM) | -0.3708782 | 2.830225 | 48.36083 | -0.1310419 | 0.9999310 | |
| Vehicle | (MSH3aso-3uM) | 0.0562575 | 2.929650 | 47.84872 | 0.0192028 | 1.0000000 | |
| Vehicle | (SCRaso-3uM) | 4.3492758 | 2.830225 | 48.36083 | 1.5367243 | 0.5442953 | |
| (MSH3aso-0.022uM) | (MSH3aso-0.26uM) | -5.3901193 | 2.986661 | 47.09787 | -1.8047309 | 0.3830559 | |
| (MSH3aso-0.022uM) | (MSH3aso-3uM) | -4.9629835 | 3.102618 | 47.32952 | -1.5996115 | 0.5050440 | |
| (MSH3aso-0.022uM) | (SCRaso-3uM) | -0.6699653 | 2.986661 | 47.09787 | -0.2243192 | 0.9994170 | |
| (MSH3aso-0.26uM) | (MSH3aso-3uM) | 0.4271357 | 3.100407 | 47.32645 | 0.1377676 | 0.9999157 | |
| (MSH3aso-0.26uM) | (SCRaso-3uM) | 4.7201540 | 2.984531 | 47.09175 | 1.5815398 | 0.5163374 | |
| (MSH3aso-3uM) | (SCRaso-3uM) | 4.2930183 | 3.100407 | 47.32645 | 1.3846626 | 0.6402912 |
Group mean at each timepoint
This performs separate ANOVA
(aov(MSH3_level ~ treatment)) and pairwise comparisons with
posthoc Tukey’s HSD tests for each time point to compare the treatment
groups. The results show pairwise comparisons of treatment group means
at each specific time point, providing insights into how treatments
differ at each stage of the experiment. The results are displayed in
separate tables for each time point.
# Vector of time points
my_timepoints <- unique(plot_data$treatment_duration_weeks)
names(my_timepoints) <- my_timepoints
# Perform ANOVA and Tukey's HSD at each time point
group_mean_compare_timepoints <- lapply(my_timepoints, function(my_timepoint,
plot_data) {
# Filter data
my_data <- plot_data %>%
dplyr::filter(treatment_duration_weeks == my_timepoint)
# ANOVA
my_aov <- aov(level ~ treatment, data = my_data)
# Posthoc tukey's HSD for pairwise comparisons
my_pwc <- TukeyHSD(my_aov)
# Tidy results
my_pwc_tidy <- tidy(my_pwc) %>%
dplyr::mutate(timepoint = my_timepoint,
p.signif = symnum(adj.p.value,
corr = FALSE,
na = FALSE,
legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", ""))) %>%
relocate(timepoint)
# Output
return(my_pwc_tidy)
},
plot_data = plot_data)
# Edit names
names(group_mean_compare_timepoints) <- paste0("week", names(group_mean_compare_timepoints))
# Display tables
for (table_name in names(group_mean_compare_timepoints)) {
print(kable(group_mean_compare_timepoints[[table_name]],
caption = paste(table_name, "pairwise comparison of treatment group means")))
}##
##
## Table: week3 pairwise comparison of treatment group means
##
## | timepoint|term |contrast | null.value| estimate| conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|---------:|----------:|---------:|-----------:|:--------|
## | 3|treatment |MSH3aso-0.022uM-Vehicle | 0| 3.75360| -110.53382| 118.04102| 0.9999662| |
## | 3|treatment |MSH3aso-0.26uM-Vehicle | 0| -80.68483| -194.97225| 33.60258| 0.2207749| |
## | 3|treatment |MSH3aso-3uM-Vehicle | 0| -94.02397| -208.31138| 20.26345| 0.1251589| |
## | 3|treatment |SCRaso-3uM-Vehicle | 0| 90.14653| -24.14088| 204.43395| 0.1482275| |
## | 3|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM | 0| -84.43843| -206.61682| 37.73995| 0.2365331| |
## | 3|treatment |MSH3aso-3uM-MSH3aso-0.022uM | 0| -97.77757| -219.95595| 24.40082| 0.1400028| |
## | 3|treatment |SCRaso-3uM-MSH3aso-0.022uM | 0| 86.39293| -35.78545| 208.57132| 0.2196206| |
## | 3|treatment |MSH3aso-3uM-MSH3aso-0.26uM | 0| -13.33913| -135.51752| 108.83925| 0.9961523| |
## | 3|treatment |SCRaso-3uM-MSH3aso-0.26uM | 0| 170.83137| 48.65298| 293.00975| 0.0061940|** |
## | 3|treatment |SCRaso-3uM-MSH3aso-3uM | 0| 184.17050| 61.99211| 306.34889| 0.0035671|** |
##
##
## Table: week9 pairwise comparison of treatment group means
##
## | timepoint|term |contrast | null.value| estimate| conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|----------:|------------:|---------:|-----------:|:--------|
## | 9|treatment |MSH3aso-0.022uM-Vehicle | 0| -49.257150| -151.4014689| 52.88717| 0.5009701| |
## | 9|treatment |MSH3aso-0.26uM-Vehicle | 0| -52.523533| -142.6063552| 37.55929| 0.3391298| |
## | 9|treatment |MSH3aso-3uM-Vehicle | 0| -98.764800| -230.6325487| 33.10295| 0.1628924| |
## | 9|treatment |SCRaso-3uM-Vehicle | 0| 59.393400| -30.6894219| 149.47622| 0.2445285| |
## | 9|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM | 0| -3.266383| -110.9359493| 104.40318| 0.9999661| |
## | 9|treatment |MSH3aso-3uM-MSH3aso-0.022uM | 0| -49.507650| -193.9615311| 94.94623| 0.7603098| |
## | 9|treatment |SCRaso-3uM-MSH3aso-0.022uM | 0| 108.650550| 0.9809841| 216.32012| 0.0479058|* |
## | 9|treatment |MSH3aso-3uM-MSH3aso-0.26uM | 0| -46.241267| -182.4336919| 89.95116| 0.7659766| |
## | 9|treatment |SCRaso-3uM-MSH3aso-0.26uM | 0| 111.916933| 15.6143459| 208.21952| 0.0235927|* |
## | 9|treatment |SCRaso-3uM-MSH3aso-3uM | 0| 158.158200| 21.9657748| 294.35063| 0.0236845|* |
##
##
## Table: week12 pairwise comparison of treatment group means
##
## | timepoint|term |contrast | null.value| estimate| conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|---------:|-----------:|---------:|-----------:|:--------|
## | 12|treatment |MSH3aso-0.022uM-Vehicle | 0| -65.72353| -112.375298| -19.07177| 0.0058676|** |
## | 12|treatment |MSH3aso-0.26uM-Vehicle | 0| -79.95663| -126.608397| -33.30487| 0.0013056|** |
## | 12|treatment |MSH3aso-3uM-Vehicle | 0| -99.32160| -145.973364| -52.66984| 0.0002043|*** |
## | 12|treatment |SCRaso-3uM-Vehicle | 0| 50.23097| 3.579203| 96.88273| 0.0333263|* |
## | 12|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM | 0| -14.23310| -64.105934| 35.63973| 0.8822168| |
## | 12|treatment |MSH3aso-3uM-MSH3aso-0.022uM | 0| -33.59807| -83.470900| 16.27477| 0.2561342| |
## | 12|treatment |SCRaso-3uM-MSH3aso-0.022uM | 0| 115.95450| 66.081666| 165.82733| 0.0000917|*** |
## | 12|treatment |MSH3aso-3uM-MSH3aso-0.26uM | 0| -19.36497| -69.237800| 30.50787| 0.7213519| |
## | 12|treatment |SCRaso-3uM-MSH3aso-0.26uM | 0| 130.18760| 80.314766| 180.06043| 0.0000309|*** |
## | 12|treatment |SCRaso-3uM-MSH3aso-3uM | 0| 149.55257| 99.679733| 199.42540| 0.0000080|*** |
##
##
## Table: week15 pairwise comparison of treatment group means
##
## | timepoint|term |contrast | null.value| estimate| conf.low| conf.high| adj.p.value|p.signif |
## |---------:|:---------|:------------------------------|----------:|---------:|-----------:|---------:|-----------:|:--------|
## | 15|treatment |MSH3aso-0.022uM-Vehicle | 0| -53.42553| -150.078839| 43.22777| 0.4144185| |
## | 15|treatment |MSH3aso-0.26uM-Vehicle | 0| -76.49517| -173.148473| 20.15814| 0.1426786| |
## | 15|treatment |MSH3aso-3uM-Vehicle | 0| -96.69235| -206.286898| 12.90220| 0.0909143|. |
## | 15|treatment |SCRaso-3uM-Vehicle | 0| 32.53630| -64.117006| 129.18961| 0.7990783| |
## | 15|treatment |MSH3aso-0.26uM-MSH3aso-0.022uM | 0| -23.06963| -126.396364| 80.25710| 0.9432480| |
## | 15|treatment |MSH3aso-3uM-MSH3aso-0.022uM | 0| -43.26682| -158.789613| 72.25598| 0.7342527| |
## | 15|treatment |SCRaso-3uM-MSH3aso-0.022uM | 0| 85.96183| -17.364897| 189.28856| 0.1169050| |
## | 15|treatment |MSH3aso-3uM-MSH3aso-0.26uM | 0| -20.19718| -135.719980| 95.32561| 0.9758222| |
## | 15|treatment |SCRaso-3uM-MSH3aso-0.26uM | 0| 109.03147| 5.704736| 212.35820| 0.0377067|* |
## | 15|treatment |SCRaso-3uM-MSH3aso-3uM | 0| 129.22865| 13.705853| 244.75145| 0.0272796|* |MSH3ko western
## [1] "MSH3ko_western_MSH2_MSH6"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!protein,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("genotype", "rep")) %>%
dplyr::mutate(level = as.numeric(level) * 100) %>%
dplyr::filter(!is.na(level))
# Relevel variables
plot_data <- plot_data %>%
dplyr::mutate(genotype = fct_relevel(genotype, intersect(genotype_order, unique(plot_data$genotype))))
# Comparisons for stat_compare_means
pairwise_comparisons <- lapply(combn(unique(plot_data$genotype), 2, simplify = FALSE), as.character)
# t-test with correction
ttest <-
plot_data %>%
group_by(protein) %>%
t_test(level ~ genotype) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# Display t-test
kbl(ttest, caption = "T-test") %>%
kable_styling(bootstrap_options = "striped")| protein | .y. | group1 | group2 | n1 | n2 | statistic | df | p | p.adj | p.adj.signif | y.position | groups |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| MSH2 | level | Unedited | MSH3ko | 12 | 9 | 1.5826844 | 18.42727 | 0.131 | 0.262 | ns | 176.8654 | Unedited, MSH3ko |
| MSH6 | level | Unedited | MSH3ko | 12 | 9 | 0.1298626 | 18.34958 | 0.898 | 0.898 | ns | 224.7324 | Unedited, MSH3ko |
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ genotype * protein, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ genotype | protein) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
group_by(protein) %>%
dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::left_join(y_max, by = join_by(protein)) %>%
group_by(protein) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
ungroup()
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
kable_styling(bootstrap_options = "striped")| contrast | group1 | group2 | protein | estimate | SE | df | t.ratio | p.value | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|---|
| Unedited - MSH3ko | Unedited | MSH3ko | MSH2 | 18.37798 | 16.42012 | 38 | 1.1192356 | 0.2700647 | 0.5401294 | ns |
| Unedited - MSH3ko | Unedited | MSH3ko | MSH6 | 2.50350 | 16.42012 | 38 | 0.1524654 | 0.8796266 | 0.8796266 | ns |
# Plot
my_plot <-
ggplot(plot_data, aes(x = genotype, y = level)) +
facet_wrap(vars(protein)) +
stat_summary(aes(fill = genotype),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = genotype),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
# geom_violin(trim = TRUE) +
# geom_boxplot(outlier.shape = NA, width = 0.2) +
# geom_jitter(aes(colour = genotype, fill = genotype),
# width = 0.2, height = 0, size = point_size, alpha = 0.75) +
# stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
# stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
# aes(label = after_stat(round(y, 1)))) +
# stat_compare_means(comparisons = pairwise_comparisons,
# label = "p.signif",
# hide.ns = hide_ns) +
# stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj") +
stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
geom_hline(yintercept = 0, colour = "darkgrey") +
geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
scale_colour_manual(values = genotype_colour) +
scale_fill_manual(values = genotype_colour) +
scale_x_discrete(labels = genotype_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Genotype",
y = "Relative protein level (%)",
colour = "Genotype",
fill = "Genotype") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
my_plot
Small pool PCR
Data from FAN1ko and unedited expansion and contraction tables has been combined into a single facetted plot.
Analysis
# Experiment IDs
experiment_ids <- setNames(names(MSH3aso_data)[6:9], names(MSH3aso_data)[6:9])
# Make results directory
dir.create(file.path(out_dir, "sppcr"),
recursive = TRUE)
# Format data
plot_data <- lapply(experiment_ids, function(experiment_id,
MSH3aso_data,
treatment_rename) {
# Extract annotation from experiment ID
experiment_id_components <- as.character(str_split(experiment_id, "_", simplify = TRUE))
my_genotype = experiment_id_components[[2]]
my_metric = experiment_id_components[[3]]
# Format
my_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!Metric,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
dplyr::mutate(level = as.numeric(level)) %>%
dplyr::filter(!is.na(level)) %>%
dplyr::mutate(genotype = my_genotype,
metric = my_metric,
genotype = ifelse(genotype == "unedited", "Unedited", genotype)) %>%
relocate(c(metric, genotype), .after = "Metric")
# Output
return(my_data)
},
MSH3aso_data = MSH3aso_data,
treatment_rename = treatment_rename)
# rbind
plot_data <- rbindlist(plot_data)
# Exclude SCR treatment from unedited genotype
plot_data <- plot_data %>%
dplyr::filter(!(genotype == "Unedited" & treatment == "SCRaso-3uM"))
# Relevel
plot_data <- plot_data %>%
dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))),
# genotype = fct_relevel(genotype, intersect(genotype_order, unique(plot_data$genotype))),
genotype = fct_relevel(genotype, c("FAN1ko", "Unedited")),
Metric = fct_relevel(Metric, c("Expansions per lane", "Contractions per lane")))
# Add to treatment_rename
my_treatment_rename <- sapply(names(treatment_rename), function(x) {
if (x == "Baseline") {
paste(treatment_rename[x], "(36d)")
} else {
paste(treatment_rename[x], "(12wks)")
}
}, USE.NAMES = TRUE)
# t-test
ttest <-
plot_data %>%
group_by(Metric, genotype) %>%
t_test(level ~ treatment) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# Display t-test
kbl(ttest, caption = "T-test") %>%
kable_styling(bootstrap_options = "striped")| Metric | genotype | .y. | group1 | group2 | n1 | n2 | statistic | df | p | p.adj | p.adj.signif | y.position | groups |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Expansions per lane | FAN1ko | level | Baseline | Vehicle | 3 | 3 | -1.6223360 | 3.145469 | 0.199 | 0.8516250 | ns | 0.9830 | Baseline, Vehicle |
| Expansions per lane | FAN1ko | level | Baseline | MSH3aso-3uM | 3 | 3 | 0.3630952 | 3.997965 | 0.735 | 0.8516250 | ns | 1.1126 | Baseline , MSH3aso-3uM |
| Expansions per lane | FAN1ko | level | Baseline | SCRaso-3uM | 3 | 3 | -1.3874776 | 3.032716 | 0.258 | 0.8516250 | ns | 1.2422 | Baseline , SCRaso-3uM |
| Expansions per lane | FAN1ko | level | Vehicle | MSH3aso-3uM | 3 | 3 | 2.1007572 | 3.188233 | 0.121 | 0.8516250 | ns | 1.3718 | Vehicle , MSH3aso-3uM |
| Expansions per lane | FAN1ko | level | Vehicle | SCRaso-3uM | 3 | 3 | 0.3786803 | 3.984831 | 0.724 | 0.8516250 | ns | 1.5014 | Vehicle , SCRaso-3uM |
| Expansions per lane | FAN1ko | level | MSH3aso-3uM | SCRaso-3uM | 3 | 3 | -1.8691854 | 3.073112 | 0.156 | 0.8516250 | ns | 1.6310 | MSH3aso-3uM, SCRaso-3uM |
| Expansions per lane | Unedited | level | Baseline | Vehicle | 3 | 3 | -0.5341094 | 3.172790 | 0.628 | 0.8516250 | ns | 1.1080 | Baseline, Vehicle |
| Expansions per lane | Unedited | level | Baseline | MSH3aso-3uM | 3 | 3 | -0.0960038 | 2.991441 | 0.930 | 0.9847059 | ns | 1.2700 | Baseline , MSH3aso-3uM |
| Expansions per lane | Unedited | level | Vehicle | MSH3aso-3uM | 3 | 3 | 0.3322949 | 3.960746 | 0.757 | 0.8516250 | ns | 1.4320 | Vehicle , MSH3aso-3uM |
| Contractions per lane | FAN1ko | level | Baseline | Vehicle | 3 | 3 | -0.6167433 | 3.079932 | 0.580 | 0.8516250 | ns | 0.6790 | Baseline, Vehicle |
| Contractions per lane | FAN1ko | level | Baseline | MSH3aso-3uM | 3 | 3 | 0.4007954 | 3.967370 | 0.709 | 0.8516250 | ns | 0.8086 | Baseline , MSH3aso-3uM |
| Contractions per lane | FAN1ko | level | Baseline | SCRaso-3uM | 3 | 3 | 0.0194998 | 3.999821 | 0.985 | 0.9850000 | ns | 0.9382 | Baseline , SCRaso-3uM |
| Contractions per lane | FAN1ko | level | Vehicle | MSH3aso-3uM | 3 | 3 | 0.8921462 | 2.922612 | 0.440 | 0.8516250 | ns | 1.0678 | Vehicle , MSH3aso-3uM |
| Contractions per lane | FAN1ko | level | Vehicle | SCRaso-3uM | 3 | 3 | 0.6289575 | 3.092149 | 0.573 | 0.8516250 | ns | 1.1974 | Vehicle , SCRaso-3uM |
| Contractions per lane | FAN1ko | level | MSH3aso-3uM | SCRaso-3uM | 3 | 3 | -0.3789710 | 3.962464 | 0.724 | 0.8516250 | ns | 1.3270 | MSH3aso-3uM, SCRaso-3uM |
| Contractions per lane | Unedited | level | Baseline | Vehicle | 3 | 3 | 0.9637981 | 2.789449 | 0.411 | 0.8516250 | ns | 0.9410 | Baseline, Vehicle |
| Contractions per lane | Unedited | level | Baseline | MSH3aso-3uM | 3 | 3 | -0.5800253 | 2.278657 | 0.614 | 0.8516250 | ns | 1.1030 | Baseline , MSH3aso-3uM |
| Contractions per lane | Unedited | level | Vehicle | MSH3aso-3uM | 3 | 3 | -0.8737264 | 2.057593 | 0.472 | 0.8516250 | ns | 1.2650 | Vehicle , MSH3aso-3uM |
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ genotype * Metric, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ genotype | Metric) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
group_by(Metric) %>%
dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::left_join(y_max, by = join_by(Metric)) %>%
group_by(Metric) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
ungroup()
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
kable_styling(bootstrap_options = "striped")| contrast | group1 | group2 | Metric | estimate | SE | df | t.ratio | p.value | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|---|
| FAN1ko - Unedited | FAN1ko | Unedited | Expansions per lane | 0.0922244 | 0.109188 | 38 | 0.8446387 | 0.4035991 | 0.8071983 | ns |
| FAN1ko - Unedited | FAN1ko | Unedited | Contractions per lane | 0.0086691 | 0.109188 | 38 | 0.0793959 | 0.9371344 | 0.9371344 | ns |
Plot together
# Plot
my_plot <-
ggplot(plot_data, aes(x = treatment, y = level)) +
facet_grid(cols = vars(Metric),
rows = vars(genotype),
labeller = labeller(genotype = genotype_rename)) +
geom_hline(yintercept = 0, colour = "darkgrey") +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
# geom_violin(trim = TRUE) +
# geom_boxplot(outlier.shape = NA, width = 0.2) +
# geom_jitter(aes(colour = treatment, fill = treatment),
# width = 0.1, height = 0, size = point_size, alpha = 0.8) +
# stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
# stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
# aes(label = after_stat(round(y, 1)))) +
stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj", tip.length = 0.01) +
scale_colour_manual(values = treatment_colour) +
scale_fill_manual(values = treatment_colour) +
scale_x_discrete(labels = my_treatment_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Treatment",
y = "Instability events per lane",
colour = "Treatment",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
my_plot
Plot genotypes separately
# Genotype vector
my_genotypes <- as.character(unique(plot_data$genotype))
names(my_genotypes) <- my_genotypes
# Plot
my_plots <- lapply(my_genotypes, function(my_genotype,
plot_data,
ttest) {
# Filter data
my_data <- plot_data %>%
dplyr::filter(genotype == my_genotype)
# Plot
my_plot <-
ggplot(my_data, aes(x = treatment, y = level)) +
facet_wrap(vars(Metric)) +
geom_hline(yintercept = 0, colour = "darkgrey") +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj", tip.length = 0.01) +
scale_colour_manual(values = treatment_colour) +
scale_fill_manual(values = treatment_colour) +
scale_x_discrete(labels = my_treatment_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Treatment",
y = "Instability events per lane",
colour = "Treatment",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
# Output
return(my_plot)
},
plot_data = plot_data,
ttest = ttest)
# Display plots
for (my_plot in my_plots) {
print(my_plot)
}
Plot genotypes and metrics separately
# Genotype vector
my_genotypes <- as.character(unique(plot_data$genotype))
names(my_genotypes) <- my_genotypes
# Metric vector
my_metrics <- as.character(unique(plot_data$Metric))
names(my_metrics) <- my_metrics
# Plot
my_plots <- lapply(my_genotypes, function(my_genotype,
plot_data,
ttest,
my_metrics) {
# Plot for each metric
my_plot_metrics <- lapply(my_metrics, function(my_metric,
my_genotype) {
# # Manual
# my_metric = my_metrics[[1]]
# Filter data
my_data <- plot_data %>%
dplyr::filter(genotype == my_genotype,
Metric == my_metric)
# Plot
my_plot <-
ggplot(my_data, aes(x = treatment, y = level)) +
geom_hline(yintercept = 0, colour = "darkgrey") +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj", tip.length = 0.01) +
scale_colour_manual(values = treatment_colour) +
scale_fill_manual(values = treatment_colour) +
scale_x_discrete(labels = my_treatment_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Treatment",
y = "Instability events per lane",
colour = "Treatment",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
# Output
return(my_plot)
},
my_genotype = my_genotype)
# Output
return(my_plot_metrics)
},
plot_data = plot_data,
ttest = ttest,
my_metrics = my_metrics)
# Flatten
my_plots <- unlist(my_plots, recursive = FALSE)
# Display plots
for (my_plot in my_plots) {
print(my_plot)
}
In vivo titration
Import data
# Experiment IDs
experiment_ids <- setNames(names(MSH3aso_data)[10:11], names(MSH3aso_data)[10:11])
# Make results directory
dir.create(file.path(out_dir, "in_vivo_titration"),
recursive = TRUE)
# Format data
plot_data <- lapply(experiment_ids, function(experiment_id,
MSH3aso_data) {
# Format
my_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!tissue,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
dplyr::mutate(level = as.numeric(level),
experiment_id = experiment_id) %>%
dplyr::filter(!is.na(level)) %>%
relocate(experiment_id)
# Output
return(my_data)
},
MSH3aso_data = MSH3aso_data)
# rbind
plot_data <- rbindlist(plot_data)
# Relevel
plot_data <- plot_data %>%
dplyr::mutate(tissue = fct_relevel(tissue, intersect(tissue_order, unique(plot_data$tissue))),
treatment = fct_relevel(treatment, intersect(treatment, unique(plot_data$treatment))))Experiments combined
Data from two experiments have been combined for stats and plotting.
# t-test
ttest <-
plot_data %>%
group_by(tissue) %>%
t_test(level ~ treatment) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# Display t-test
kbl(ttest, caption = "T-test") %>%
kable_styling(bootstrap_options = "striped")| tissue | .y. | group1 | group2 | n1 | n2 | statistic | df | p | p.adj | p.adj.signif | y.position | groups |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Cortex | level | Vehicle | 3ug | 8 | 8 | 1.5415216 | 9.132123 | 1.57e-01 | 0.1860741 | ns | 119.3469 | Vehicle, 3ug |
| Cortex | level | Vehicle | 10ug | 8 | 8 | 2.5865314 | 8.569527 | 3.00e-02 | 0.0417391 |
|
130.4767 | Vehicle, 10ug |
| Cortex | level | Vehicle | 30ug | 8 | 8 | 13.0916812 | 9.896674 | 1.00e-07 | 0.0000023 | *** | 141.6066 | Vehicle, 30ug |
| Cortex | level | Vehicle | SCRaso | 8 | 4 | -0.8948261 | 5.464449 | 4.09e-01 | 0.4362667 | ns | 152.7365 | Vehicle, SCRaso |
| Cortex | level | 3ug | 10ug | 8 | 8 | 1.0006693 | 13.661257 | 3.34e-01 | 0.3685517 | ns | 163.8663 | 3ug , 10ug |
| Cortex | level | 3ug | 30ug | 8 | 8 | 8.0797152 | 13.639093 | 1.50e-06 | 0.0000088 | *** | 174.9962 | 3ug , 30ug |
| Cortex | level | 3ug | SCRaso | 8 | 4 | -1.9654142 | 9.987477 | 7.80e-02 | 0.0998400 | ns | 186.1261 | 3ug , SCRaso |
| Cortex | level | 10ug | 30ug | 8 | 8 | 6.2577056 | 12.757057 | 3.20e-05 | 0.0000853 | *** | 197.2559 | 10ug, 30ug |
| Cortex | level | 10ug | SCRaso | 8 | 4 | -2.9044349 | 9.731249 | 1.60e-02 | 0.0232727 |
|
208.3858 | 10ug , SCRaso |
| Cortex | level | 30ug | SCRaso | 8 | 4 | -12.2273936 | 9.858507 | 3.00e-07 | 0.0000030 | *** | 219.5157 | 30ug , SCRaso |
| Striatum | level | Vehicle | 3ug | 4 | 4 | 0.5318878 | 4.900537 | 6.18e-01 | 0.6180000 | ns | 120.4869 | Vehicle, 3ug |
| Striatum | level | Vehicle | 10ug | 4 | 4 | 2.2498355 | 3.606236 | 9.50e-02 | 0.1169231 | ns | 132.5071 | Vehicle, 10ug |
| Striatum | level | Vehicle | 30ug | 4 | 5 | 7.1738608 | 4.011652 | 2.00e-03 | 0.0035556 | ** | 144.5274 | Vehicle, 30ug |
| Striatum | level | 3ug | 10ug | 4 | 4 | 2.5708958 | 4.585510 | 5.40e-02 | 0.0720000 | ns | 156.5476 | 3ug , 10ug |
| Striatum | level | 3ug | 30ug | 4 | 5 | 9.8404489 | 5.570527 | 9.93e-05 | 0.0002444 | *** | 168.5679 | 3ug , 30ug |
| Striatum | level | 10ug | 30ug | 4 | 5 | 10.3652603 | 6.925037 | 1.82e-05 | 0.0000529 | *** | 180.5882 | 10ug, 30ug |
| Brainstem | level | Vehicle | 3ug | 4 | 4 | 6.2762166 | 5.062854 | 1.00e-03 | 0.0018824 | ** | 113.6829 | Vehicle, 3ug |
| Brainstem | level | Vehicle | 10ug | 4 | 4 | 17.6634697 | 5.767722 | 3.00e-06 | 0.0000122 | *** | 125.7031 | Vehicle, 10ug |
| Brainstem | level | Vehicle | 30ug | 4 | 5 | 32.4960641 | 6.975689 | 0.00e+00 | 0.0000002 | *** | 137.7234 | Vehicle, 30ug |
| Brainstem | level | 3ug | 10ug | 4 | 4 | 8.0782464 | 5.643734 | 2.62e-04 | 0.0005989 | *** | 149.7436 | 3ug , 10ug |
| Brainstem | level | 3ug | 30ug | 4 | 5 | 18.8693737 | 5.569056 | 2.90e-06 | 0.0000122 | *** | 161.7639 | 3ug , 30ug |
| Brainstem | level | 10ug | 30ug | 4 | 5 | 12.2559102 | 6.544504 | 9.40e-06 | 0.0000299 | *** | 173.7842 | 10ug, 30ug |
| Spinal cord | level | Vehicle | 3ug | 8 | 8 | 3.1598645 | 9.857447 | 1.00e-02 | 0.0152381 |
|
117.3129 | Vehicle, 3ug |
| Spinal cord | level | Vehicle | 10ug | 8 | 8 | 5.6296674 | 7.974333 | 4.99e-04 | 0.0009980 | *** | 128.4427 | Vehicle, 10ug |
| Spinal cord | level | Vehicle | 30ug | 8 | 8 | 12.3280829 | 8.050690 | 1.70e-06 | 0.0000088 | *** | 139.5726 | Vehicle, 30ug |
| Spinal cord | level | Vehicle | SCRaso | 8 | 4 | 1.4272505 | 3.660822 | 2.33e-01 | 0.2662857 | ns | 150.7025 | Vehicle, SCRaso |
| Spinal cord | level | 3ug | 10ug | 8 | 8 | 3.3243802 | 11.143146 | 7.00e-03 | 0.0112000 |
|
161.8323 | 3ug , 10ug |
| Spinal cord | level | 3ug | 30ug | 8 | 8 | 9.2092605 | 11.401101 | 1.30e-06 | 0.0000088 | *** | 172.9622 | 3ug , 30ug |
| Spinal cord | level | 3ug | SCRaso | 8 | 4 | -0.7888956 | 6.115341 | 4.60e-01 | 0.4748387 | ns | 184.0921 | 3ug , SCRaso |
| Spinal cord | level | 10ug | 30ug | 8 | 8 | 4.6709975 | 13.979700 | 3.62e-04 | 0.0007723 | *** | 195.2219 | 10ug, 30ug |
| Spinal cord | level | 10ug | SCRaso | 8 | 4 | -3.5884061 | 9.592936 | 5.00e-03 | 0.0084211 | ** | 206.3518 | 10ug , SCRaso |
| Spinal cord | level | 30ug | SCRaso | 8 | 4 | -8.8324037 | 9.429196 | 7.30e-06 | 0.0000260 | *** | 217.4817 | 30ug , SCRaso |
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * tissue, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
group_by(tissue) %>%
dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::left_join(y_max, by = join_by(tissue)) %>%
group_by(tissue) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
ungroup()
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
kable_styling(bootstrap_options = "striped")| contrast | group1 | group2 | tissue | estimate | SE | df | t.ratio | p.value | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|---|
| Vehicle - 3ug | Vehicle | 3ug | Cortex | 6.688529 | 5.038409 | 88 | 1.3275080 | 0.6749243 | 0.7999102 | ns |
| Vehicle - 10ug | Vehicle | 10ug | Cortex | 12.910300 | 5.038409 | 88 | 2.5623763 | 0.0864607 | 0.1317497 | ns |
| Vehicle - 30ug | Vehicle | 30ug | Cortex | 49.453178 | 5.038409 | 88 | 9.8152365 | 0.0000000 | 0.0000000 | *** |
| Vehicle - SCRaso | Vehicle | SCRaso | Cortex | -2.683435 | 6.170766 | 88 | -0.4348625 | 0.9924383 | 0.9924383 | ns |
| 3ug - 10ug | 3ug | 10ug | Cortex | 6.221771 | 5.038409 | 88 | 1.2348682 | 0.7310900 | 0.8355315 | ns |
| 3ug - 30ug | 3ug | 30ug | Cortex | 42.764649 | 5.038409 | 88 | 8.4877285 | 0.0000000 | 0.0000000 | *** |
| 3ug - SCRaso | 3ug | SCRaso | Cortex | -9.371964 | 6.170766 | 88 | -1.5187683 | 0.5532082 | 0.6808716 | ns |
| 10ug - 30ug | 10ug | 30ug | Cortex | 36.542878 | 5.038409 | 88 | 7.2528603 | 0.0000000 | 0.0000000 | *** |
| 10ug - SCRaso | 10ug | SCRaso | Cortex | -15.593735 | 6.170766 | 88 | -2.5270340 | 0.0938909 | 0.1365686 | ns |
| 30ug - SCRaso | 30ug | SCRaso | Cortex | -52.136613 | 6.170766 | 88 | -8.4489696 | 0.0000000 | 0.0000000 | *** |
| Vehicle - 3ug | Vehicle | 3ug | Striatum | 3.659846 | 7.125387 | 88 | 0.5136347 | 0.9556332 | 0.9864601 | ns |
| Vehicle - 10ug | Vehicle | 10ug | Striatum | 13.950394 | 7.125387 | 88 | 1.9578438 | 0.2119660 | 0.2826214 | ns |
| Vehicle - 30ug | Vehicle | 30ug | Striatum | 45.805947 | 6.759735 | 88 | 6.7762930 | 0.0000000 | 0.0000000 | *** |
| Vehicle - SCRaso | Vehicle | SCRaso | Striatum | NA | NA | NA | NA | NA | NA | |
| 3ug - 10ug | 3ug | 10ug | Striatum | 10.290548 | 7.125387 | 88 | 1.4442091 | 0.4755368 | 0.6086871 | ns |
| 3ug - 30ug | 3ug | 30ug | Striatum | 42.146101 | 6.759735 | 88 | 6.2348745 | 0.0000001 | 0.0000002 | *** |
| 3ug - SCRaso | 3ug | SCRaso | Striatum | NA | NA | NA | NA | NA | NA | |
| 10ug - 30ug | 10ug | 30ug | Striatum | 31.855553 | 6.759735 | 88 | 4.7125445 | 0.0000532 | 0.0001064 | *** |
| 10ug - SCRaso | 10ug | SCRaso | Striatum | NA | NA | NA | NA | NA | NA | |
| 30ug - SCRaso | 30ug | SCRaso | Striatum | NA | NA | NA | NA | NA | NA | |
| Vehicle - 3ug | Vehicle | 3ug | Brainstem | 18.896578 | 7.125387 | 88 | 2.6520074 | 0.0459968 | 0.0735949 | ns |
| Vehicle - 10ug | Vehicle | 10ug | Brainstem | 44.900223 | 7.125387 | 88 | 6.3014438 | 0.0000001 | 0.0000002 | *** |
| Vehicle - 30ug | Vehicle | 30ug | Brainstem | 77.182277 | 6.759735 | 88 | 11.4179438 | 0.0000000 | 0.0000000 | *** |
| Vehicle - SCRaso | Vehicle | SCRaso | Brainstem | NA | NA | NA | NA | NA | NA | |
| 3ug - 10ug | 3ug | 10ug | Brainstem | 26.003645 | 7.125387 | 88 | 3.6494364 | 0.0024713 | 0.0041622 | ** |
| 3ug - 30ug | 3ug | 30ug | Brainstem | 58.285700 | 6.759735 | 88 | 8.6224826 | 0.0000000 | 0.0000000 | *** |
| 3ug - SCRaso | 3ug | SCRaso | Brainstem | NA | NA | NA | NA | NA | NA | |
| 10ug - 30ug | 10ug | 30ug | Brainstem | 32.282054 | 6.759735 | 88 | 4.7756388 | 0.0000416 | 0.0000887 | *** |
| 10ug - SCRaso | 10ug | SCRaso | Brainstem | NA | NA | NA | NA | NA | NA | |
| 30ug - SCRaso | 30ug | SCRaso | Brainstem | NA | NA | NA | NA | NA | NA | |
| Vehicle - 3ug | Vehicle | 3ug | Spinal cord | 11.934323 | 5.038409 | 88 | 2.3686688 | 0.1337486 | 0.1860850 | ns |
| Vehicle - 10ug | Vehicle | 10ug | Spinal cord | 34.876693 | 5.038409 | 88 | 6.9221637 | 0.0000000 | 0.0000000 | *** |
| Vehicle - 30ug | Vehicle | 30ug | Spinal cord | 73.707582 | 5.038409 | 88 | 14.6291378 | 0.0000000 | 0.0000000 | *** |
| Vehicle - SCRaso | Vehicle | SCRaso | Spinal cord | 7.259351 | 6.170766 | 88 | 1.1764101 | 0.7647539 | 0.8438664 | ns |
| 3ug - 10ug | 3ug | 10ug | Spinal cord | 22.942370 | 5.038409 | 88 | 4.5534949 | 0.0001613 | 0.0003037 | *** |
| 3ug - 30ug | 3ug | 30ug | Spinal cord | 61.773259 | 5.038409 | 88 | 12.2604689 | 0.0000000 | 0.0000000 | *** |
| 3ug - SCRaso | 3ug | SCRaso | Spinal cord | -4.674972 | 6.170766 | 88 | -0.7576000 | 0.9419017 | 0.9864601 | ns |
| 10ug - 30ug | 10ug | 30ug | Spinal cord | 38.830889 | 5.038409 | 88 | 7.7069741 | 0.0000000 | 0.0000000 | *** |
| 10ug - SCRaso | 10ug | SCRaso | Spinal cord | -27.617342 | 6.170766 | 88 | -4.4755130 | 0.0002166 | 0.0003851 | *** |
| 30ug - SCRaso | 30ug | SCRaso | Spinal cord | -66.448231 | 6.170766 | 88 | -10.7682309 | 0.0000000 | 0.0000000 | *** |
# Shape rename
my_shape_rename <- setNames(c("1", "2"), unique(plot_data$experiment_id))
# Plot
my_plot <-
ggplot(plot_data, aes(x = treatment, y = level)) +
facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment, shape = experiment_id),
width = 0.2, height = 0, size = point_size, colour = "black") +
# geom_violin(trim = TRUE) +
# geom_boxplot(outlier.shape = NA, width = 0.2) +
# geom_jitter(aes(colour = treatment, fill = treatment, shape = experiment_id),
# width = 0.1, height = 0, size = point_size, alpha = 0.8) +
# stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
# stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
# aes(label = after_stat(round(y, 1)))) +
geom_hline(yintercept = 0, colour = "darkgrey") +
geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
# stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = TRUE, tip.length = 0.01) +
scale_x_discrete(labels = treatment_ionis_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
scale_colour_manual(values = treatment_ionis_colour,
labels = treatment_ionis_rename) +
scale_fill_manual(values = treatment_ionis_colour,
labels = treatment_ionis_rename) +
# scale_shape_discrete(labels = my_shape_rename) +
scale_shape_manual(values = c(21, 24),
labels = my_shape_rename) +
labs(x = "Treatment",
y = "Relative expression (%)",
colour = "Treatment",
fill = "Treatment",
shape = "Experiment") +
theme_minimal() +
theme(text = element_text(size = element_text_size))
my_plot
Experiment separated
Data from the two experiments kept separate for stats and plotting.
# Vector of experiments
my_experiment_names <- unique(plot_data$experiment_id)
names(my_experiment_names) <- my_experiment_names
# Plot
my_plots <- lapply(my_experiment_names, function(my_experiment_name,
plot_data,
padj_method,
hide_ns,
treatment_ionis_rename,
treatment_ionis_colour,
element_text_size) {
# Filter data
my_data <- plot_data %>%
dplyr::filter(experiment_id == my_experiment_name)
# t-test
ttest <-
my_data %>%
group_by(tissue) %>%
t_test(level ~ treatment) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * tissue, data = my_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- my_data %>%
group_by(tissue) %>%
dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::left_join(y_max, by = join_by(tissue)) %>%
group_by(tissue) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
ungroup()
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Plot
my_plot <-
ggplot(my_data, aes(x = treatment, y = level)) +
facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
geom_hline(yintercept = 0, colour = "darkgrey") +
geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
# stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = hide_ns, tip.length = 0.01) +
scale_x_discrete(labels = treatment_ionis_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
scale_colour_manual(values = treatment_ionis_colour) +
scale_fill_manual(values = treatment_ionis_colour) +
labs(x = "Treatment",
y = "Relative expression (%)",
colour = "Treatment",
fill = "Treatment") +
theme_minimal() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
# Output
out <- mget(c("ttest", "my_pwc", "my_plot"))
return(out)
},
plot_data = plot_data,
padj_method = padj_method,
hide_ns = hide_ns,
treatment_ionis_rename = treatment_ionis_rename,
treatment_ionis_colour = treatment_ionis_colour,
element_text_size = element_text_size)
# Display ttests
my_ttests <- lapply(my_experiment_names, function(my_experiment_name,
my_plots) {
my_plots[[my_experiment_name]]$ttest
},
my_plots = my_plots)
for (my_experiment_name in my_experiment_names) {
my_ttest <- my_ttests[[my_experiment_name]]
print(kable(my_ttest, caption = paste(my_experiment_name, "t-test")))
}##
##
## Table: invivo_titration t-test
##
## |tissue |.y. |group1 |group2 | n1| n2| statistic| df| p| p.adj|p.adj.signif | y.position|groups |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|--------:|---------:|:------------|----------:|:-------------|
## |Cortex |level |Vehicle |3ug | 4| 4| 1.8542000| 3.934168| 1.39e-01| 0.1516364|ns | 116.4359|Vehicle, 3ug |
## |Cortex |level |Vehicle |10ug | 4| 4| 1.9604778| 3.502115| 1.32e-01| 0.1508571|ns | 128.4561|Vehicle, 10ug |
## |Cortex |level |Vehicle |30ug | 4| 5| 12.0544775| 6.802275| 7.70e-06| 0.0000449|*** | 140.4764|Vehicle, 30ug |
## |Cortex |level |3ug |10ug | 4| 4| 0.4778913| 5.476986| 6.51e-01| 0.6510000|ns | 152.4966|3ug , 10ug |
## |Cortex |level |3ug |30ug | 4| 5| 5.4406421| 4.771674| 3.00e-03| 0.0055385|** | 164.5169|3ug , 30ug |
## |Cortex |level |10ug |30ug | 4| 5| 3.6419864| 3.964163| 2.20e-02| 0.0368000|* | 176.5372|10ug, 30ug |
## |Striatum |level |Vehicle |3ug | 4| 4| 0.5318878| 4.900537| 6.18e-01| 0.6448696|ns | 120.4869|Vehicle, 3ug |
## |Striatum |level |Vehicle |10ug | 4| 4| 2.2498355| 3.606236| 9.50e-02| 0.1200000|ns | 132.5071|Vehicle, 10ug |
## |Striatum |level |Vehicle |30ug | 4| 5| 7.1738608| 4.011652| 2.00e-03| 0.0040000|** | 144.5274|Vehicle, 30ug |
## |Striatum |level |3ug |10ug | 4| 4| 2.5708958| 4.585510| 5.40e-02| 0.0720000|ns | 156.5476|3ug , 10ug |
## |Striatum |level |3ug |30ug | 4| 5| 9.8404489| 5.570527| 9.93e-05| 0.0002648|*** | 168.5679|3ug , 30ug |
## |Striatum |level |10ug |30ug | 4| 5| 10.3652603| 6.925037| 1.82e-05| 0.0000624|*** | 180.5882|10ug, 30ug |
## |Brainstem |level |Vehicle |3ug | 4| 4| 6.2762166| 5.062854| 1.00e-03| 0.0021818|** | 113.6829|Vehicle, 3ug |
## |Brainstem |level |Vehicle |10ug | 4| 4| 17.6634697| 5.767722| 3.00e-06| 0.0000244|*** | 125.7031|Vehicle, 10ug |
## |Brainstem |level |Vehicle |30ug | 4| 5| 32.4960641| 6.975689| 0.00e+00| 0.0000002|*** | 137.7234|Vehicle, 30ug |
## |Brainstem |level |3ug |10ug | 4| 4| 8.0782464| 5.643734| 2.62e-04| 0.0006288|*** | 149.7436|3ug , 10ug |
## |Brainstem |level |3ug |30ug | 4| 5| 18.8693737| 5.569056| 2.90e-06| 0.0000244|*** | 161.7639|3ug , 30ug |
## |Brainstem |level |10ug |30ug | 4| 5| 12.2559102| 6.544504| 9.40e-06| 0.0000449|*** | 173.7842|10ug, 30ug |
## |Spinal cord |level |Vehicle |3ug | 4| 4| 1.8137090| 5.196380| 1.27e-01| 0.1508571|ns | 114.0809|Vehicle, 3ug |
## |Spinal cord |level |Vehicle |10ug | 4| 4| 4.1886198| 3.144929| 2.30e-02| 0.0368000|* | 126.1011|Vehicle, 10ug |
## |Spinal cord |level |Vehicle |30ug | 4| 5| 15.5784053| 4.806628| 2.68e-05| 0.0000804|*** | 138.1214|Vehicle, 30ug |
## |Spinal cord |level |3ug |10ug | 4| 4| 3.6277004| 3.331969| 3.00e-02| 0.0423529|* | 150.1416|3ug , 10ug |
## |Spinal cord |level |3ug |30ug | 4| 5| 13.7512712| 5.706043| 1.37e-05| 0.0000548|*** | 162.1619|3ug , 30ug |
## |Spinal cord |level |10ug |30ug | 4| 5| 3.2870097| 4.377864| 2.60e-02| 0.0390000|* | 174.1822|10ug, 30ug |
##
##
## Table: invivo_titration_SCR t-test
##
## |tissue |.y. |group1 |group2 | n1| n2| statistic| df| p| p.adj|p.adj.signif | y.position|groups |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|-----:|---------:|:------------|----------:|:---------------|
## |Cortex |level |Vehicle |3ug | 4| 4| 0.1049662| 5.594147| 0.920| 0.9200000|ns | 115.1301|Vehicle, 3ug |
## |Cortex |level |Vehicle |10ug | 4| 4| 2.2245897| 5.918406| 0.068| 0.1133333|ns | 121.5746|Vehicle, 10ug |
## |Cortex |level |Vehicle |30ug | 4| 3| 6.7692358| 2.581470| 0.010| 0.0500000|* | 128.0191|Vehicle, 30ug |
## |Cortex |level |Vehicle |SCRaso | 4| 4| -0.8455827| 5.916147| 0.431| 0.4600000|ns | 134.4637|Vehicle, SCRaso |
## |Cortex |level |3ug |10ug | 4| 4| 1.8321285| 5.855849| 0.118| 0.1653333|ns | 140.9082|3ug , 10ug |
## |Cortex |level |3ug |30ug | 4| 3| 6.4166254| 3.000136| 0.008| 0.0500000|* | 147.3527|3ug , 30ug |
## |Cortex |level |3ug |SCRaso | 4| 4| -0.8348014| 5.858741| 0.437| 0.4600000|ns | 153.7973|3ug , SCRaso |
## |Cortex |level |10ug |30ug | 4| 3| 5.4902697| 2.734738| 0.015| 0.0500000|* | 160.2418|10ug, 30ug |
## |Cortex |level |10ug |SCRaso | 4| 4| -2.9027231| 5.999983| 0.027| 0.0700000|ns | 166.6863|10ug , SCRaso |
## |Cortex |level |30ug |SCRaso | 3| 4| -7.1014215| 2.737148| 0.008| 0.0500000|* | 173.1309|30ug , SCRaso |
## |Spinal cord |level |Vehicle |3ug | 4| 4| 3.5857849| 5.441242| 0.014| 0.0500000|* | 113.0961|Vehicle, 3ug |
## |Spinal cord |level |Vehicle |10ug | 4| 4| 4.4059859| 4.737829| 0.008| 0.0500000|* | 119.5406|Vehicle, 10ug |
## |Spinal cord |level |Vehicle |30ug | 4| 3| 4.8340610| 2.218758| 0.033| 0.0700000|ns | 125.9851|Vehicle, 30ug |
## |Spinal cord |level |Vehicle |SCRaso | 4| 4| 1.2633565| 5.058685| 0.262| 0.3082353|ns | 132.4297|Vehicle, SCRaso |
## |Spinal cord |level |3ug |10ug | 4| 4| 1.2454440| 5.687328| 0.262| 0.3082353|ns | 138.8742|3ug , 10ug |
## |Spinal cord |level |3ug |30ug | 4| 3| 3.3573428| 2.426350| 0.060| 0.1090909|ns | 145.3187|3ug , 30ug |
## |Spinal cord |level |3ug |SCRaso | 4| 4| -1.7908931| 5.902224| 0.124| 0.1653333|ns | 151.7633|3ug , SCRaso |
## |Spinal cord |level |10ug |30ug | 4| 3| 2.6587774| 2.686702| 0.086| 0.1323077|ns | 158.2078|10ug, 30ug |
## |Spinal cord |level |10ug |SCRaso | 4| 4| -2.7716089| 5.929482| 0.033| 0.0700000|ns | 164.6523|10ug , SCRaso |
## |Spinal cord |level |30ug |SCRaso | 3| 4| -4.1376927| 2.552372| 0.035| 0.0700000|ns | 171.0969|30ug , SCRaso |# Display ANOVA posthoc
my_pwcs <- lapply(my_experiment_names, function(my_experiment_name,
my_plots) {
my_plots[[my_experiment_name]]$my_pwc
},
my_plots = my_plots)
for (my_experiment_name in my_experiment_names) {
my_pwc <- my_pwcs[[my_experiment_name]]
print(kable(my_pwc, caption = paste(my_experiment_name, "ANOVA with corrected posthoc pairwise comparisons")))
}##
##
## Table: invivo_titration ANOVA with corrected posthoc pairwise comparisons
##
## |contrast |group1 |group2 |tissue | estimate| SE| df| t.ratio| p.value| p.adj|p.adj.signif |
## |:--------------|:-------|:------|:-----------|---------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 3ug |Vehicle |3ug |Cortex | 12.985657| 6.591570| 52| 1.9700402| 0.2125051| 0.2550061|ns |
## |Vehicle - 10ug |Vehicle |10ug |Cortex | 18.273183| 6.591570| 52| 2.7722052| 0.0374764| 0.0499685|* |
## |Vehicle - 30ug |Vehicle |30ug |Cortex | 53.393389| 6.253312| 52| 8.5384178| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Cortex | 5.287527| 6.591570| 52| 0.8021650| 0.8530927| 0.9016176|ns |
## |3ug - 30ug |3ug |30ug |Cortex | 40.407732| 6.253312| 52| 6.4618131| 0.0000002| 0.0000005|*** |
## |10ug - 30ug |10ug |30ug |Cortex | 35.120206| 6.253312| 52| 5.6162569| 0.0000045| 0.0000084|*** |
## |Vehicle - 3ug |Vehicle |3ug |Striatum | 3.659846| 6.591570| 52| 0.5552313| 0.9447105| 0.9447105|ns |
## |Vehicle - 10ug |Vehicle |10ug |Striatum | 13.950394| 6.591570| 52| 2.1163994| 0.1614185| 0.2038971|ns |
## |Vehicle - 30ug |Vehicle |30ug |Striatum | 45.805947| 6.253312| 52| 7.3250699| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Striatum | 10.290548| 6.591570| 52| 1.5611681| 0.4094088| 0.4678958|ns |
## |3ug - 30ug |3ug |30ug |Striatum | 42.146101| 6.253312| 52| 6.7398047| 0.0000001| 0.0000002|*** |
## |10ug - 30ug |10ug |30ug |Striatum | 31.855553| 6.253312| 52| 5.0941891| 0.0000288| 0.0000461|*** |
## |Vehicle - 3ug |Vehicle |3ug |Brainstem | 18.896578| 6.591570| 52| 2.8667797| 0.0295351| 0.0416967|* |
## |Vehicle - 10ug |Vehicle |10ug |Brainstem | 44.900223| 6.591570| 52| 6.8117651| 0.0000001| 0.0000002|*** |
## |Vehicle - 30ug |Vehicle |30ug |Brainstem | 77.182277| 6.253312| 52| 12.3426240| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Brainstem | 26.003645| 6.591570| 52| 3.9449854| 0.0013358| 0.0020037|** |
## |3ug - 30ug |3ug |30ug |Brainstem | 58.285700| 6.253312| 52| 9.3207729| 0.0000000| 0.0000000|*** |
## |10ug - 30ug |10ug |30ug |Brainstem | 32.282054| 6.253312| 52| 5.1623931| 0.0000227| 0.0000389|*** |
## |Vehicle - 3ug |Vehicle |3ug |Spinal cord | 5.128342| 6.591570| 52| 0.7780154| 0.8640502| 0.9016176|ns |
## |Vehicle - 10ug |Vehicle |10ug |Spinal cord | 42.464438| 6.591570| 52| 6.4422348| 0.0000002| 0.0000005|*** |
## |Vehicle - 30ug |Vehicle |30ug |Spinal cord | 79.016220| 6.253312| 52| 12.6358995| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Spinal cord | 37.336096| 6.591570| 52| 5.6642194| 0.0000038| 0.0000076|*** |
## |3ug - 30ug |3ug |30ug |Spinal cord | 73.887878| 6.253312| 52| 11.8157992| 0.0000000| 0.0000000|*** |
## |10ug - 30ug |10ug |30ug |Spinal cord | 36.551782| 6.253312| 52| 5.8451877| 0.0000020| 0.0000043|*** |
##
##
## Table: invivo_titration_SCR ANOVA with corrected posthoc pairwise comparisons
##
## |contrast |group1 |group2 |tissue | estimate| SE| df| t.ratio| p.value| p.adj|p.adj.signif |
## |:----------------|:-------|:------|:-----------|-----------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 3ug |Vehicle |3ug |Cortex | 0.3914004| 6.702412| 28| 0.0583970| 0.9999972| 0.9999972|ns |
## |Vehicle - 10ug |Vehicle |10ug |Cortex | 7.5474167| 6.702412| 28| 1.1260747| 0.7915231| 0.9661809|ns |
## |Vehicle - 30ug |Vehicle |30ug |Cortex | 43.0115304| 7.239432| 28| 5.9412850| 0.0000200| 0.0000771|*** |
## |Vehicle - SCRaso |Vehicle |SCRaso |Cortex | -2.8714926| 6.702412| 28| -0.4284268| 0.9925764| 0.9999972|ns |
## |3ug - 10ug |3ug |10ug |Cortex | 7.1560163| 6.702412| 28| 1.0676778| 0.8212538| 0.9661809|ns |
## |3ug - 30ug |3ug |30ug |Cortex | 42.6201300| 7.239432| 28| 5.8872199| 0.0000231| 0.0000771|*** |
## |3ug - SCRaso |3ug |SCRaso |Cortex | -3.2628930| 6.702412| 28| -0.4868237| 0.9879628| 0.9999972|ns |
## |10ug - 30ug |10ug |30ug |Cortex | 35.4641137| 7.239432| 28| 4.8987424| 0.0003293| 0.0008232|*** |
## |10ug - SCRaso |10ug |SCRaso |Cortex | -10.4189093| 6.702412| 28| -1.5545015| 0.5374236| 0.8268055|ns |
## |30ug - SCRaso |30ug |SCRaso |Cortex | -45.8830230| 7.239432| 28| -6.3379311| 0.0000070| 0.0000349|*** |
## |Vehicle - 3ug |Vehicle |3ug |Spinal cord | 18.7403032| 6.702412| 28| 2.7960536| 0.0646660| 0.1175745|ns |
## |Vehicle - 10ug |Vehicle |10ug |Spinal cord | 27.2889479| 6.702412| 28| 4.0715116| 0.0029549| 0.0065663|** |
## |Vehicle - 30ug |Vehicle |30ug |Spinal cord | 64.8874611| 7.239432| 28| 8.9630593| 0.0000000| 0.0000002|*** |
## |Vehicle - SCRaso |Vehicle |SCRaso |Spinal cord | 7.2179372| 6.702412| 28| 1.0769164| 0.8166817| 0.9661809|ns |
## |3ug - 10ug |3ug |10ug |Spinal cord | 8.5486447| 6.702412| 28| 1.2754580| 0.7078962| 0.9661809|ns |
## |3ug - 30ug |3ug |30ug |Spinal cord | 46.1471579| 7.239432| 28| 6.3744167| 0.0000063| 0.0000349|*** |
## |3ug - SCRaso |3ug |SCRaso |Spinal cord | -11.5223660| 6.702412| 28| -1.7191372| 0.4390542| 0.7317570|ns |
## |10ug - 30ug |10ug |30ug |Spinal cord | 37.5985132| 7.239432| 28| 5.1935720| 0.0001490| 0.0004257|*** |
## |10ug - SCRaso |10ug |SCRaso |Spinal cord | -20.0710107| 6.702412| 28| -2.9945952| 0.0416897| 0.0833794|ns |
## |30ug - SCRaso |30ug |SCRaso |Spinal cord | -57.6695239| 7.239432| 28| -7.9660285| 0.0000001| 0.0000011|*** |# Display plots
my_bars <- lapply(my_experiment_names, function(my_experiment_name,
my_plots) {
my_plots[[my_experiment_name]]$my_plot
},
my_plots = my_plots)
for (my_experiment_name in my_experiment_names) {
print(my_experiment_name)
print(my_bars[[my_experiment_name]])
}## [1] "invivo_titration"
## [1] "invivo_titration_SCR"
# Export plots
# for (my_experiment_name in my_experiment_names) {
# for (output_format in output_formats) {
# ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration", "_", my_experiment_name, ".", output_format)),
# plot = my_bars[[my_experiment_name]],
# height = 6, width = 12)
# }
# }
# Export experiment 1 plot
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-experiment1", ".", output_format)),
plot = my_bars[[my_experiment_names[[1]]]],
height = 6, width = 12)
}
# Export experiment 2 plot
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-experiment2", ".", output_format)),
plot = my_bars[[my_experiment_names[[2]]]],
height = 6, width = 9)
}Combined experiments with SCR separate
# List of plots
my_treatment_groups <- list(MSH3aso = c("Vehicle", "3ug", "10ug", "30ug"),
SCRaso = c("Vehicle", "SCRaso"))
# Plot
my_plots <- lapply(my_treatment_groups, function(my_treatment_group,
plot_data,
padj_method,
hide_ns,
treatment_ionis_rename,
treatment_ionis_colour,
element_text_size) {
# Filter data
my_data <- plot_data %>%
dplyr::filter(treatment %in% my_treatment_group)
# Filter to keep only tissues where all treatments are present
my_data <- my_data %>%
group_by(tissue) %>%
filter(all(my_treatment_group %in% treatment)) %>%
ungroup()
# t-test
ttest <-
my_data %>%
group_by(tissue) %>%
t_test(level ~ treatment) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * tissue, data = my_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- my_data %>%
group_by(tissue) %>%
dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::left_join(y_max, by = join_by(tissue)) %>%
group_by(tissue) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
ungroup()
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Plot
my_plot <-
ggplot(my_data, aes(x = treatment, y = level)) +
facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
geom_hline(yintercept = 0, colour = "darkgrey") +
geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
# stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = hide_ns, tip.length = 0.01) +
scale_x_discrete(labels = treatment_ionis_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
scale_colour_manual(values = treatment_ionis_colour) +
scale_fill_manual(values = treatment_ionis_colour) +
labs(x = "Treatment",
y = "Relative expression (%)",
colour = "Treatment",
fill = "Treatment") +
theme_minimal() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
# Output
out <- mget(c("ttest", "my_pwc", "my_plot"))
return(out)
},
plot_data = plot_data,
padj_method = padj_method,
hide_ns = hide_ns,
treatment_ionis_rename = treatment_ionis_rename,
treatment_ionis_colour = treatment_ionis_colour,
element_text_size = element_text_size)
# Vector of treatment group names
my_treatment_group_names <- setNames(names(my_treatment_groups), names(my_treatment_groups))
# Display ttests
my_ttests <- lapply(my_treatment_group_names, function(my_treatment_group_name,
my_plots) {
my_plots[[my_treatment_group_name]]$ttest
},
my_plots = my_plots)
for (my_treatment_group_name in my_treatment_group_names) {
my_ttest <- my_ttests[[my_treatment_group_name]]
print(kable(my_ttest, caption = paste(my_treatment_group_name, "t-test")))
}##
##
## Table: MSH3aso t-test
##
## |tissue |.y. |group1 |group2 | n1| n2| statistic| df| p| p.adj|p.adj.signif | y.position|groups |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|---------:|--------:|---------:|:------------|----------:|:-------------|
## |Cortex |level |Vehicle |3ug | 8| 8| 1.5415216| 9.132123| 1.57e-01| 0.1712727|ns | 116.4359|Vehicle, 3ug |
## |Cortex |level |Vehicle |10ug | 8| 8| 2.5865314| 8.569527| 3.00e-02| 0.0378947|* | 128.4561|Vehicle, 10ug |
## |Cortex |level |Vehicle |30ug | 8| 8| 13.0916812| 9.896674| 1.00e-07| 0.0000017|*** | 140.4764|Vehicle, 30ug |
## |Cortex |level |3ug |10ug | 8| 8| 1.0006693| 13.661257| 3.34e-01| 0.3485217|ns | 152.4966|3ug , 10ug |
## |Cortex |level |3ug |30ug | 8| 8| 8.0797152| 13.639093| 1.50e-06| 0.0000079|*** | 164.5169|3ug , 30ug |
## |Cortex |level |10ug |30ug | 8| 8| 6.2577056| 12.757057| 3.20e-05| 0.0000768|*** | 176.5372|10ug, 30ug |
## |Striatum |level |Vehicle |3ug | 4| 4| 0.5318878| 4.900537| 6.18e-01| 0.6180000|ns | 120.4869|Vehicle, 3ug |
## |Striatum |level |Vehicle |10ug | 4| 4| 2.2498355| 3.606236| 9.50e-02| 0.1085714|ns | 132.5071|Vehicle, 10ug |
## |Striatum |level |Vehicle |30ug | 4| 5| 7.1738608| 4.011652| 2.00e-03| 0.0030000|** | 144.5274|Vehicle, 30ug |
## |Striatum |level |3ug |10ug | 4| 4| 2.5708958| 4.585510| 5.40e-02| 0.0648000|ns | 156.5476|3ug , 10ug |
## |Striatum |level |3ug |30ug | 4| 5| 9.8404489| 5.570527| 9.93e-05| 0.0002167|*** | 168.5679|3ug , 30ug |
## |Striatum |level |10ug |30ug | 4| 5| 10.3652603| 6.925037| 1.82e-05| 0.0000485|*** | 180.5882|10ug, 30ug |
## |Brainstem |level |Vehicle |3ug | 4| 4| 6.2762166| 5.062854| 1.00e-03| 0.0016000|** | 113.6829|Vehicle, 3ug |
## |Brainstem |level |Vehicle |10ug | 4| 4| 17.6634697| 5.767722| 3.00e-06| 0.0000105|*** | 125.7031|Vehicle, 10ug |
## |Brainstem |level |Vehicle |30ug | 4| 5| 32.4960641| 6.975689| 0.00e+00| 0.0000002|*** | 137.7234|Vehicle, 30ug |
## |Brainstem |level |3ug |10ug | 4| 4| 8.0782464| 5.643734| 2.62e-04| 0.0005240|*** | 149.7436|3ug , 10ug |
## |Brainstem |level |3ug |30ug | 4| 5| 18.8693737| 5.569056| 2.90e-06| 0.0000105|*** | 161.7639|3ug , 30ug |
## |Brainstem |level |10ug |30ug | 4| 5| 12.2559102| 6.544504| 9.40e-06| 0.0000281|*** | 173.7842|10ug, 30ug |
## |Spinal cord |level |Vehicle |3ug | 8| 8| 3.1598645| 9.857447| 1.00e-02| 0.0133333|* | 117.3129|Vehicle, 3ug |
## |Spinal cord |level |Vehicle |10ug | 8| 8| 5.6296674| 7.974333| 4.99e-04| 0.0008554|*** | 129.3331|Vehicle, 10ug |
## |Spinal cord |level |Vehicle |30ug | 8| 8| 12.3280829| 8.050690| 1.70e-06| 0.0000079|*** | 141.3534|Vehicle, 30ug |
## |Spinal cord |level |3ug |10ug | 8| 8| 3.3243802| 11.143146| 7.00e-03| 0.0098824|** | 153.3736|3ug , 10ug |
## |Spinal cord |level |3ug |30ug | 8| 8| 9.2092605| 11.401101| 1.30e-06| 0.0000079|*** | 165.3939|3ug , 30ug |
## |Spinal cord |level |10ug |30ug | 8| 8| 4.6709975| 13.979700| 3.62e-04| 0.0006683|*** | 177.4142|10ug, 30ug |
##
##
## Table: SCRaso t-test
##
## |tissue |.y. |group1 |group2 | n1| n2| statistic| df| p| p.adj|p.adj.signif | y.position|groups |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|-----:|-----:|:------------|----------:|:---------------|
## |Cortex |level |Vehicle |SCRaso | 8| 4| -0.8948261| 5.464449| 0.409| 0.409|ns | 110.1899|Vehicle, SCRaso |
## |Spinal cord |level |Vehicle |SCRaso | 8| 4| 1.4272505| 3.660822| 0.233| 0.409|ns | 108.1559|Vehicle, SCRaso |# Display ANOVA posthoc
my_pwcs <- lapply(my_treatment_group_names, function(my_treatment_group_name,
my_plots) {
my_plots[[my_treatment_group_name]]$my_pwc
},
my_plots = my_plots)
for (my_treatment_group_name in my_treatment_group_names) {
my_pwc <- my_pwcs[[my_treatment_group_name]]
print(kable(my_pwc, caption = paste(my_treatment_group_name, "ANOVA with corrected posthoc pairwise comparisons")))
}##
##
## Table: MSH3aso ANOVA with corrected posthoc pairwise comparisons
##
## |contrast |group1 |group2 |tissue | estimate| SE| df| t.ratio| p.value| p.adj|p.adj.signif |
## |:--------------|:-------|:------|:-----------|---------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 3ug |Vehicle |3ug |Cortex | 6.688529| 5.113932| 82| 1.3079033| 0.5605969| 0.6115603|ns |
## |Vehicle - 10ug |Vehicle |10ug |Cortex | 12.910300| 5.113932| 82| 2.5245348| 0.0635297| 0.0847063|ns |
## |Vehicle - 30ug |Vehicle |30ug |Cortex | 49.453178| 5.113932| 82| 9.6702841| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Cortex | 6.221771| 5.113932| 82| 1.2166316| 0.6181520| 0.6450282|ns |
## |3ug - 30ug |3ug |30ug |Cortex | 42.764649| 5.113932| 82| 8.3623808| 0.0000000| 0.0000000|*** |
## |10ug - 30ug |10ug |30ug |Cortex | 36.542878| 5.113932| 82| 7.1457493| 0.0000000| 0.0000000|*** |
## |Vehicle - 3ug |Vehicle |3ug |Striatum | 3.659846| 7.232192| 82| 0.5060493| 0.9574291| 0.9574291|ns |
## |Vehicle - 10ug |Vehicle |10ug |Striatum | 13.950394| 7.232192| 82| 1.9289301| 0.2240923| 0.2689108|ns |
## |Vehicle - 30ug |Vehicle |30ug |Striatum | 45.805947| 6.861060| 82| 6.6762200| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Striatum | 10.290548| 7.232192| 82| 1.4228808| 0.4888740| 0.5587132|ns |
## |3ug - 30ug |3ug |30ug |Striatum | 42.146101| 6.861060| 82| 6.1427972| 0.0000002| 0.0000003|*** |
## |10ug - 30ug |10ug |30ug |Striatum | 31.855553| 6.861060| 82| 4.6429491| 0.0000752| 0.0001289|*** |
## |Vehicle - 3ug |Vehicle |3ug |Brainstem | 18.896578| 7.232192| 82| 2.6128422| 0.0512146| 0.0723029|ns |
## |Vehicle - 10ug |Vehicle |10ug |Brainstem | 44.900223| 7.232192| 82| 6.2083834| 0.0000001| 0.0000003|*** |
## |Vehicle - 30ug |Vehicle |30ug |Brainstem | 77.182277| 6.861060| 82| 11.2493224| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Brainstem | 26.003645| 7.232192| 82| 3.5955412| 0.0030423| 0.0045635|** |
## |3ug - 30ug |3ug |30ug |Brainstem | 58.285700| 6.861060| 82| 8.4951449| 0.0000000| 0.0000000|*** |
## |10ug - 30ug |10ug |30ug |Brainstem | 32.282054| 6.861060| 82| 4.7051117| 0.0000593| 0.0001095|*** |
## |Vehicle - 3ug |Vehicle |3ug |Spinal cord | 11.934323| 5.113932| 82| 2.3336881| 0.0987803| 0.1247751|ns |
## |Vehicle - 10ug |Vehicle |10ug |Spinal cord | 34.876693| 5.113932| 82| 6.8199365| 0.0000000| 0.0000000|*** |
## |Vehicle - 30ug |Vehicle |30ug |Spinal cord | 73.707582| 5.113932| 82| 14.4130931| 0.0000000| 0.0000000|*** |
## |3ug - 10ug |3ug |10ug |Spinal cord | 22.942370| 5.113932| 82| 4.4862484| 0.0001357| 0.0002171|*** |
## |3ug - 30ug |3ug |30ug |Spinal cord | 61.773259| 5.113932| 82| 12.0794051| 0.0000000| 0.0000000|*** |
## |10ug - 30ug |10ug |30ug |Spinal cord | 38.830889| 5.113932| 82| 7.5931567| 0.0000000| 0.0000000|*** |
##
##
## Table: SCRaso ANOVA with corrected posthoc pairwise comparisons
##
## |contrast |group1 |group2 |tissue | estimate| SE| df| t.ratio| p.value| p.adj|p.adj.signif |
## |:----------------|:-------|:------|:-----------|---------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - SCRaso |Vehicle |SCRaso |Cortex | -2.683435| 3.465514| 20| -0.7743251| 0.4477962| 0.4477962|ns |
## |Vehicle - SCRaso |Vehicle |SCRaso |Spinal cord | 7.259351| 3.465514| 20| 2.0947398| 0.0491315| 0.0982631|ns |# Display plots
my_bars <- lapply(my_treatment_group_names, function(my_treatment_group_name,
my_plots) {
my_plots[[my_treatment_group_name]]$my_plot
},
my_plots = my_plots)
for (my_treatment_group_name in my_treatment_group_names) {
print(my_treatment_group_name)
print(my_bars[[my_treatment_group_name]])
}## [1] "MSH3aso"
## [1] "SCRaso"
# Export plots
# for (my_experiment_name in my_experiment_names) {
# for (output_format in output_formats) {
# ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration", "_", my_experiment_name, ".", output_format)),
# plot = my_bars[[my_experiment_name]],
# height = 6, width = 12)
# }
# }
# Export experiment 1 plot
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-combined-", my_treatment_group_names[[1]], ".", output_format)),
plot = my_bars[[my_treatment_group_names[[1]]]],
height = 6, width = 12)
}
# Export experiment 2 plot
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, "in_vivo_titration", paste0("in_vivo_titration-combined-", my_treatment_group_names[[2]], ".", output_format)),
plot = my_bars[[my_treatment_group_names[[2]]]],
height = 5, width = 7)
}Custom
Experiment 2 with only Vehicle, SCR ASO and the highest dose of MSH3 ASO.
# Filter data
my_data <- plot_data %>%
dplyr::filter(experiment_id == "invivo_titration_SCR",
treatment %in% c("Vehicle", "SCRaso", "30ug"))
# t-test
ttest <-
my_data %>%
group_by(tissue) %>%
t_test(level ~ treatment) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# Display ttest
print(kable(ttest, caption = "t-test"))##
##
## Table: t-test
##
## |tissue |.y. |group1 |group2 | n1| n2| statistic| df| p| p.adj|p.adj.signif | y.position|groups |
## |:-----------|:-----|:-------|:------|--:|--:|----------:|--------:|-----:|------:|:------------|----------:|:---------------|
## |Cortex |level |Vehicle |30ug | 4| 3| 6.7692358| 2.581470| 0.010| 0.0300|* | 115.1301|Vehicle, 30ug |
## |Cortex |level |Vehicle |SCRaso | 4| 4| -0.8455827| 5.916147| 0.431| 0.4310|ns | 123.8302|Vehicle, SCRaso |
## |Cortex |level |30ug |SCRaso | 3| 4| -7.1014215| 2.737148| 0.008| 0.0300|* | 132.5303|30ug , SCRaso |
## |Spinal cord |level |Vehicle |30ug | 4| 3| 4.8340610| 2.218758| 0.033| 0.0525|ns | 113.0961|Vehicle, 30ug |
## |Spinal cord |level |Vehicle |SCRaso | 4| 4| 1.2633565| 5.058685| 0.262| 0.3144|ns | 121.7962|Vehicle, SCRaso |
## |Spinal cord |level |30ug |SCRaso | 3| 4| -4.1376927| 2.552372| 0.035| 0.0525|ns | 130.4963|30ug , SCRaso |# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * tissue, data = my_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | tissue) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- my_data %>%
group_by(tissue) %>%
dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::left_join(y_max, by = join_by(tissue)) %>%
group_by(tissue) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.2 * unique(ymax), length.out = n())) %>%
ungroup()
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Display ANOVA
print(kable(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons"))##
##
## Table: ANOVA with corrected posthoc pairwise comparisons
##
## |contrast |group1 |group2 |tissue | estimate| SE| df| t.ratio| p.value| p.adj|p.adj.signif |
## |:----------------|:-------|:------|:-----------|----------:|--------:|--:|----------:|---------:|---------:|:------------|
## |Vehicle - 30ug |Vehicle |30ug |Cortex | 43.011530| 8.026822| 16| 5.3584755| 0.0001792| 0.0002688|*** |
## |Vehicle - SCRaso |Vehicle |SCRaso |Cortex | -2.871493| 7.431393| 16| -0.3864003| 0.9213944| 0.9213944|ns |
## |30ug - SCRaso |30ug |SCRaso |Cortex | -45.883023| 8.026822| 16| -5.7162127| 0.0000896| 0.0001792|*** |
## |Vehicle - 30ug |Vehicle |30ug |Spinal cord | 64.887461| 8.026822| 16| 8.0838294| 0.0000014| 0.0000083|*** |
## |Vehicle - SCRaso |Vehicle |SCRaso |Spinal cord | 7.217937| 7.431393| 16| 0.9712764| 0.6047812| 0.7257374|ns |
## |30ug - SCRaso |30ug |SCRaso |Spinal cord | -57.669524| 8.026822| 16| -7.1846021| 0.0000062| 0.0000186|*** |# Plot
my_plot <-
ggplot(my_data, aes(x = treatment, y = level)) +
facet_wrap(vars(tissue), nrow = 1, scales = "free_x") +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
geom_hline(yintercept = 0, colour = "darkgrey") +
geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
# stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj.signif")
stat_pvalue_manual(my_pwc_plot, label = "p.adj.signif", hide.ns = hide_ns, tip.length = 0.01) +
scale_x_discrete(labels = treatment_ionis_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
scale_colour_manual(values = treatment_ionis_colour) +
scale_fill_manual(values = treatment_ionis_colour) +
labs(x = "Treatment",
y = "Relative expression (%)",
colour = "Treatment",
fill = "Treatment") +
theme_minimal() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
# Display plot
print(my_plot)
ASO screen
Fit models
## [1] "aso_screen"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(!log_conc_nM,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("treatment", "metric")) %>%
pivot_wider(names_from = "metric",
values_from = "level") %>%
dplyr::mutate(log_conc_nM = as.numeric(log_conc_nM),
level = as.numeric(level)) %>%
dplyr::mutate(conc_nM = exp(log_conc_nM)) %>%
relocate(conc_nM) %>%
dplyr::filter(!is.na(level))
# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
models = my_models,
data = plot_data,
xvar = "conc_nM",
xlabel = "MSH3 ASO concentration (nM)",
yvar = "level",
ylabel = "Relative expression (%)",
group_var = "treatment",
group_label = "Treatment",
log_increment = log_increment,
plot_se = FALSE)
# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
my_fit$my_plot
})
wrap_plots(my_plots) +
plot_annotation(title = "Regression models",
theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))
# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
r2 = unlist(my_r2)) %>%
tibble::rownames_to_column("model_number") %>%
arrange(-r2)
# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
kable_styling(bootstrap_options = "striped")| model_number | model_name | r2 |
|---|---|---|
| 4 | poly3 | 1.0000000 |
| 3 | poly2 | 0.9906616 |
| 2 | log | 0.9801475 |
| 1 | linear | 0.9579879 |
| 5 | exponential | 0.8681231 |
| 6 | exponential | 0.8681231 |
Plot
# Plot
my_plot <-
ggplot(plot_data, aes(x = conc_nM, y = level, colour = treatment, fill = treatment)) +
geom_point(size = point_size) +
geom_errorbar(aes(ymin = level - sd,
ymax = level + sd),
width = 0.25) +
stat_smooth(method = lm,
formula = y ~ log(x),
fullrange = FALSE, alpha = 0.2, se = FALSE) +
scale_x_continuous(breaks = pretty_breaks()) +
scale_y_continuous(trans = "log2",
breaks = c(0.1, 1, 2, 5, 10, 25, 50, 100)) +
labs(x = "MSH3 ASO dose (nM)",
y = "Relative expression (%)",
colour = "ASO",
fill = "ASO") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
plot = my_plot)
}## Saving 7 x 5 in image
Viability
## [1] "viability"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
dplyr::select(-metric) %>%
pivot_longer(!treatment_duration_weeks,
names_to = "label",
values_to = "level") %>%
separate_wider_delim(label, "_", names = c("treatment", "rep")) %>%
dplyr::mutate(treatment_duration_weeks = as.numeric(treatment_duration_weeks),
level = as.numeric(level)) %>%
dplyr::filter(!is.na(level)) %>%
dplyr::mutate(treatment_duration_weeks = case_when(treatment_duration_weeks == 9 ~ "9 weeks",
treatment_duration_weeks == 15 ~ "15 weeks",
TRUE ~ NA))
# Relevel variables
plot_data <- plot_data %>%
dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))),
treatment_duration_weeks = fct_relevel(treatment_duration_weeks, c("9 weeks", "15 weeks")))
# t-test with correction
ttest <-
plot_data %>%
group_by(treatment_duration_weeks) %>%
t_test(level ~ treatment) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# Display t-test
kbl(ttest, caption = "T-test") %>%
kable_styling(bootstrap_options = "striped")| treatment_duration_weeks | .y. | group1 | group2 | n1 | n2 | statistic | df | p | p.adj | p.adj.signif | y.position | groups |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 9 weeks | level | Vehicle | MSH3aso-0.022uM | 3 | 3 | -2.2792048 | 2.000000 | 0.150 | 0.4860000 | ns | 185.7200 | Vehicle , MSH3aso-0.022uM |
| 9 weeks | level | Vehicle | MSH3aso-0.26uM | 3 | 3 | -2.4600632 | 2.000000 | 0.133 | 0.4860000 | ns | 202.5629 | Vehicle , MSH3aso-0.26uM |
| 9 weeks | level | Vehicle | MSH3aso-3uM | 3 | 3 | -2.1686631 | 2.000000 | 0.162 | 0.4860000 | ns | 219.4057 | Vehicle , MSH3aso-3uM |
| 9 weeks | level | Vehicle | SCRaso-3uM | 3 | 3 | -2.2294431 | 2.000000 | 0.156 | 0.4860000 | ns | 236.2486 | Vehicle , SCRaso-3uM |
| 9 weeks | level | Vehicle | H2O2-500uM | 3 | 3 | 4.3821316 | 2.000000 | 0.048 | 0.2400000 | ns | 253.0914 | Vehicle , H2O2-500uM |
| 9 weeks | level | MSH3aso-0.022uM | MSH3aso-0.26uM | 3 | 3 | -0.4005964 | 3.897133 | 0.710 | 0.9800000 | ns | 269.9343 | MSH3aso-0.022uM, MSH3aso-0.26uM |
| 9 weeks | level | MSH3aso-0.022uM | MSH3aso-3uM | 3 | 3 | -0.2154569 | 3.866679 | 0.840 | 0.9800000 | ns | 286.7771 | MSH3aso-0.022uM, MSH3aso-3uM |
| 9 weeks | level | MSH3aso-0.022uM | SCRaso-3uM | 3 | 3 | -0.5169365 | 3.583734 | 0.635 | 0.9800000 | ns | 303.6200 | MSH3aso-0.022uM, SCRaso-3uM |
| 9 weeks | level | MSH3aso-0.022uM | H2O2-500uM | 3 | 3 | 4.7208869 | 3.999167 | 0.009 | 0.1125000 | ns | 320.4629 | MSH3aso-0.022uM, H2O2-500uM |
| 9 weeks | level | MSH3aso-0.26uM | MSH3aso-3uM | 3 | 3 | 0.1668522 | 3.997708 | 0.876 | 0.9800000 | ns | 337.3057 | MSH3aso-0.26uM, MSH3aso-3uM |
| 9 weeks | level | MSH3aso-0.26uM | SCRaso-3uM | 3 | 3 | -0.1521099 | 3.862324 | 0.887 | 0.9800000 | ns | 354.1486 | MSH3aso-0.26uM, SCRaso-3uM |
| 9 weeks | level | MSH3aso-0.26uM | H2O2-500uM | 3 | 3 | 4.7236434 | 3.913832 | 0.010 | 0.1125000 | ns | 370.9914 | MSH3aso-0.26uM, H2O2-500uM |
| 9 weeks | level | MSH3aso-3uM | SCRaso-3uM | 3 | 3 | -0.3012376 | 3.893224 | 0.779 | 0.9800000 | ns | 387.8343 | MSH3aso-3uM, SCRaso-3uM |
| 9 weeks | level | MSH3aso-3uM | H2O2-500uM | 3 | 3 | 4.4800005 | 3.885439 | 0.012 | 0.1125000 | ns | 404.6771 | MSH3aso-3uM, H2O2-500uM |
| 9 weeks | level | SCRaso-3uM | H2O2-500uM | 3 | 3 | 4.3568231 | 3.611482 | 0.015 | 0.1125000 | ns | 421.5200 | SCRaso-3uM, H2O2-500uM |
| 15 weeks | level | Vehicle | MSH3aso-0.022uM | 2 | 2 | -1.1058058 | 1.000277 | 0.468 | 0.9360000 | ns | 205.7200 | Vehicle , MSH3aso-0.022uM |
| 15 weeks | level | Vehicle | MSH3aso-0.26uM | 2 | 2 | -0.8782982 | 1.000365 | 0.541 | 0.9547059 | ns | 222.5629 | Vehicle , MSH3aso-0.26uM |
| 15 weeks | level | Vehicle | MSH3aso-3uM | 2 | 2 | -0.6412665 | 1.000236 | 0.637 | 0.9800000 | ns | 239.4057 | Vehicle , MSH3aso-3uM |
| 15 weeks | level | Vehicle | SCRaso-3uM | 2 | 2 | -0.9314195 | 1.000375 | 0.523 | 0.9547059 | ns | 256.2486 | Vehicle , SCRaso-3uM |
| 15 weeks | level | Vehicle | H2O2-500uM | 2 | 2 | 14.6315099 | 1.055513 | 0.038 | 0.2280000 | ns | 273.0914 | Vehicle , H2O2-500uM |
| 15 weeks | level | MSH3aso-0.022uM | MSH3aso-0.26uM | 2 | 2 | 0.2573233 | 1.962780 | 0.821 | 0.9800000 | ns | 289.9343 | MSH3aso-0.022uM, MSH3aso-0.26uM |
| 15 weeks | level | MSH3aso-0.022uM | MSH3aso-3uM | 2 | 2 | 0.2794283 | 1.987604 | 0.806 | 0.9800000 | ns | 306.7771 | MSH3aso-0.022uM, MSH3aso-3uM |
| 15 weeks | level | MSH3aso-0.022uM | SCRaso-3uM | 2 | 2 | 0.2320504 | 1.955399 | 0.839 | 0.9800000 | ns | 323.6200 | MSH3aso-0.022uM, SCRaso-3uM |
| 15 weeks | level | MSH3aso-0.022uM | H2O2-500uM | 2 | 2 | 2.1475974 | 1.009965 | 0.275 | 0.6530769 | ns | 340.4629 | MSH3aso-0.022uM, H2O2-500uM |
| 15 weeks | level | MSH3aso-0.26uM | MSH3aso-3uM | 2 | 2 | 0.0508183 | 1.912147 | 0.964 | 0.9800000 | ns | 357.3057 | MSH3aso-0.26uM, MSH3aso-3uM |
| 15 weeks | level | MSH3aso-0.26uM | SCRaso-3uM | 2 | 2 | -0.0288608 | 1.999630 | 0.980 | 0.9800000 | ns | 374.1486 | MSH3aso-0.26uM, SCRaso-3uM |
| 15 weeks | level | MSH3aso-0.26uM | H2O2-500uM | 2 | 2 | 2.0742740 | 1.013148 | 0.283 | 0.6530769 | ns | 390.9914 | MSH3aso-0.26uM, H2O2-500uM |
| 15 weeks | level | MSH3aso-3uM | SCRaso-3uM | 2 | 2 | -0.0766325 | 1.901755 | 0.946 | 0.9800000 | ns | 407.8343 | MSH3aso-3uM, SCRaso-3uM |
| 15 weeks | level | MSH3aso-3uM | H2O2-500uM | 2 | 2 | 1.6052854 | 1.008506 | 0.353 | 0.7564286 | ns | 424.6771 | MSH3aso-3uM, H2O2-500uM |
| 15 weeks | level | SCRaso-3uM | H2O2-500uM | 2 | 2 | 2.1434571 | 1.013510 | 0.275 | 0.6530769 | ns | 441.5200 | SCRaso-3uM, H2O2-500uM |
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(level ~ treatment * treatment_duration_weeks, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment | treatment_duration_weeks) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- plot_data %>%
group_by(treatment_duration_weeks) %>%
dplyr::summarise(ymax = max(level, na.rm = TRUE))
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::left_join(y_max, by = join_by(treatment_duration_weeks)) %>%
group_by(treatment_duration_weeks) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n())) %>%
ungroup()
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
kable_styling(bootstrap_options = "striped")| contrast | group1 | group2 | treatment_duration_weeks | estimate | SE | df | t.ratio | p.value | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|---|
| Vehicle - (MSH3aso-0.022uM) | Vehicle | (MSH3aso-0.022uM) | 9 weeks | -27.000000 | 26.58425 | 18 | -1.0156390 | 0.9066705 | 1.0000000 | ns |
| Vehicle - (MSH3aso-0.26uM) | Vehicle | (MSH3aso-0.26uM) | 9 weeks | -34.333333 | 26.58425 | 18 | -1.2914916 | 0.7857650 | 1.0000000 | ns |
| Vehicle - (MSH3aso-3uM) | Vehicle | (MSH3aso-3uM) | 9 weeks | -31.000000 | 26.58425 | 18 | -1.1661041 | 0.8468018 | 1.0000000 | ns |
| Vehicle - (SCRaso-3uM) | Vehicle | (SCRaso-3uM) | 9 weeks | -37.666667 | 26.58425 | 18 | -1.4168792 | 0.7170104 | 1.0000000 | ns |
| Vehicle - (H2O2-500uM) | Vehicle | (H2O2-500uM) | 9 weeks | 52.666667 | 26.58425 | 18 | 1.9811231 | 0.3896084 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (MSH3aso-0.26uM) | (MSH3aso-0.022uM) | (MSH3aso-0.26uM) | 9 weeks | -7.333333 | 26.58425 | 18 | -0.2758526 | 0.9997426 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (MSH3aso-3uM) | (MSH3aso-0.022uM) | (MSH3aso-3uM) | 9 weeks | -4.000000 | 26.58425 | 18 | -0.1504650 | 0.9999871 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (SCRaso-3uM) | (MSH3aso-0.022uM) | (SCRaso-3uM) | 9 weeks | -10.666667 | 26.58425 | 18 | -0.4012401 | 0.9984202 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (H2O2-500uM) | (MSH3aso-0.022uM) | (H2O2-500uM) | 9 weeks | 79.666667 | 26.58425 | 18 | 2.9967621 | 0.0711926 | 0.5339445 | ns |
| (MSH3aso-0.26uM) - (MSH3aso-3uM) | (MSH3aso-0.26uM) | (MSH3aso-3uM) | 9 weeks | 3.333333 | 26.58425 | 18 | 0.1253875 | 0.9999948 | 1.0000000 | ns |
| (MSH3aso-0.26uM) - (SCRaso-3uM) | (MSH3aso-0.26uM) | (SCRaso-3uM) | 9 weeks | -3.333333 | 26.58425 | 18 | -0.1253875 | 0.9999948 | 1.0000000 | ns |
| (MSH3aso-0.26uM) - (H2O2-500uM) | (MSH3aso-0.26uM) | (H2O2-500uM) | 9 weeks | 87.000000 | 26.58425 | 18 | 3.2726147 | 0.0414274 | 0.5313155 | ns |
| (MSH3aso-3uM) - (SCRaso-3uM) | (MSH3aso-3uM) | (SCRaso-3uM) | 9 weeks | -6.666667 | 26.58425 | 18 | -0.2507751 | 0.9998387 | 1.0000000 | ns |
| (MSH3aso-3uM) - (H2O2-500uM) | (MSH3aso-3uM) | (H2O2-500uM) | 9 weeks | 83.666667 | 26.58425 | 18 | 3.1472272 | 0.0531315 | 0.5313155 | ns |
| (SCRaso-3uM) - (H2O2-500uM) | (SCRaso-3uM) | (H2O2-500uM) | 9 weeks | 90.333333 | 26.58425 | 18 | 3.3980023 | 0.0321778 | 0.5313155 | ns |
| Vehicle - (MSH3aso-0.022uM) | Vehicle | (MSH3aso-0.022uM) | 15 weeks | -47.000000 | 32.55892 | 18 | -1.4435368 | 0.7016682 | 1.0000000 | ns |
| Vehicle - (MSH3aso-0.26uM) | Vehicle | (MSH3aso-0.26uM) | 15 weeks | -32.500000 | 32.55892 | 18 | -0.9981903 | 0.9125629 | 1.0000000 | ns |
| Vehicle - (MSH3aso-3uM) | Vehicle | (MSH3aso-3uM) | 15 weeks | -29.500000 | 32.55892 | 18 | -0.9060497 | 0.9399425 | 1.0000000 | ns |
| Vehicle - (SCRaso-3uM) | Vehicle | (SCRaso-3uM) | 15 weeks | -34.000000 | 32.55892 | 18 | -1.0442607 | 0.8965202 | 1.0000000 | ns |
| Vehicle - (H2O2-500uM) | Vehicle | (H2O2-500uM) | 15 weeks | 44.500000 | 32.55892 | 18 | 1.3667529 | 0.7452448 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (MSH3aso-0.26uM) | (MSH3aso-0.022uM) | (MSH3aso-0.26uM) | 15 weeks | 14.500000 | 32.55892 | 18 | 0.4453465 | 0.9974065 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (MSH3aso-3uM) | (MSH3aso-0.022uM) | (MSH3aso-3uM) | 15 weeks | 17.500000 | 32.55892 | 18 | 0.5374871 | 0.9937575 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (SCRaso-3uM) | (MSH3aso-0.022uM) | (SCRaso-3uM) | 15 weeks | 13.000000 | 32.55892 | 18 | 0.3992761 | 0.9984568 | 1.0000000 | ns |
| (MSH3aso-0.022uM) - (H2O2-500uM) | (MSH3aso-0.022uM) | (H2O2-500uM) | 15 weeks | 91.500000 | 32.55892 | 18 | 2.8102897 | 0.1012131 | 0.6072789 | ns |
| (MSH3aso-0.26uM) - (MSH3aso-3uM) | (MSH3aso-0.26uM) | (MSH3aso-3uM) | 15 weeks | 3.000000 | 32.55892 | 18 | 0.0921406 | 0.9999989 | 1.0000000 | ns |
| (MSH3aso-0.26uM) - (SCRaso-3uM) | (MSH3aso-0.26uM) | (SCRaso-3uM) | 15 weeks | -1.500000 | 32.55892 | 18 | -0.0460703 | 1.0000000 | 1.0000000 | ns |
| (MSH3aso-0.26uM) - (H2O2-500uM) | (MSH3aso-0.26uM) | (H2O2-500uM) | 15 weeks | 77.000000 | 32.55892 | 18 | 2.3649432 | 0.2200472 | 0.9430596 | ns |
| (MSH3aso-3uM) - (SCRaso-3uM) | (MSH3aso-3uM) | (SCRaso-3uM) | 15 weeks | -4.500000 | 32.55892 | 18 | -0.1382110 | 0.9999915 | 1.0000000 | ns |
| (MSH3aso-3uM) - (H2O2-500uM) | (MSH3aso-3uM) | (H2O2-500uM) | 15 weeks | 74.000000 | 32.55892 | 18 | 2.2728026 | 0.2548095 | 0.9555356 | ns |
| (SCRaso-3uM) - (H2O2-500uM) | (SCRaso-3uM) | (H2O2-500uM) | 15 weeks | 78.500000 | 32.55892 | 18 | 2.4110136 | 0.2040760 | 0.9430596 | ns |
# Plot
my_plot <-
ggplot(plot_data, aes(x = treatment, y = level)) +
facet_wrap(vars(treatment_duration_weeks)) +
stat_summary(aes(fill = treatment),
# fun = mean, colour = "black", geom = "bar", alpha = 0.8) +
fun = mean, colour = "black", geom = "bar") +
stat_summary(fun.data = mean_se,
geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
width = 0.2, height = 0, size = point_size, colour = "black", shape = 21) +
# geom_violin(trim = TRUE) +
# geom_boxplot(outlier.shape = NA, width = 0.2) +
# geom_jitter(aes(colour = treatment, fill = treatment),
# width = 0.2, height = 0, size = point_size, alpha = 0.75) +
# stat_summary(fun = mean, geom = "point", shape = 18, size = point_size + 2, show.legend = FALSE) +
# stat_summary(fun = mean, geom = "text", show.legend = FALSE, vjust = -1, size = point_size + 2,
# aes(label = after_stat(round(y, 1)))) +
# stat_compare_means(comparisons = pairwise_comparisons,
# label = "p.signif",
# hide.ns = hide_ns) +
# stat_pvalue_manual(ttest, hide.ns = hide_ns, label = "p.adj") +
stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
geom_hline(yintercept = 0, colour = "darkgrey") +
geom_hline(yintercept = 100, colour = "darkgrey", linetype = "dashed") +
scale_colour_manual(values = treatment_colour) +
scale_fill_manual(values = treatment_colour) +
scale_x_discrete(labels = \(x) scales::label_parse()(treatment_rename_subscript[x]),
guide = guide_axis(angle = 55)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Treatment",
y = "Cell viability (%)",
colour = "Treatment",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
my_plot
Mihai MSH3 normalised
## [1] "mihai_MSH3_normalised"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
pivot_longer(everything(),
names_to = "treatment",
values_to = "MSH3_counts_normalised") %>%
dplyr::mutate(MSH3_counts_normalised = as.numeric(MSH3_counts_normalised)) %>%
dplyr::filter(!is.na(MSH3_counts_normalised))
# Relevel variables
plot_data <- plot_data %>%
dplyr::mutate(treatment = fct_relevel(treatment, intersect(treatment_order, unique(plot_data$treatment))))
# t-test with correction
ttest <-
plot_data %>%
t_test(MSH3_counts_normalised ~ treatment) %>%
adjust_pvalue(method = padj_method) %>%
add_significance(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")) %>%
add_y_position() %>%
ungroup()
# Display t-test
kbl(ttest, caption = "T-test") %>%
kable_styling(bootstrap_options = "striped")| .y. | group1 | group2 | n1 | n2 | statistic | df | p | p.adj | p.adj.signif | y.position | groups |
|---|---|---|---|---|---|---|---|---|---|---|---|
| MSH3_counts_normalised | Vehicle | SCRaso | 9 | 8 | 2.015853 | 14.428636 | 6.3e-02 | 0.0630000 | ns | 722.0400 | Vehicle, SCRaso |
| MSH3_counts_normalised | Vehicle | MSH3aso_1161149 | 9 | 9 | 22.278108 | 8.165391 | 0.0e+00 | 0.0000001 | *** | 804.3067 | Vehicle , MSH3aso_1161149 |
| MSH3_counts_normalised | Vehicle | MSH3aso_1161173 | 9 | 9 | 20.646141 | 8.425659 | 0.0e+00 | 0.0000001 | *** | 886.5733 | Vehicle , MSH3aso_1161173 |
| MSH3_counts_normalised | Vehicle | MSH3aso_1161329 | 9 | 9 | 13.425819 | 15.976367 | 0.0e+00 | 0.0000000 | *** | 968.8400 | Vehicle , MSH3aso_1161329 |
| MSH3_counts_normalised | SCRaso | MSH3aso_1161149 | 8 | 9 | 16.864602 | 7.111053 | 5.0e-07 | 0.0000011 | *** | 1051.1067 | SCRaso , MSH3aso_1161149 |
| MSH3_counts_normalised | SCRaso | MSH3aso_1161173 | 8 | 9 | 15.483349 | 7.286119 | 8.0e-07 | 0.0000013 | *** | 1133.3733 | SCRaso , MSH3aso_1161173 |
| MSH3_counts_normalised | SCRaso | MSH3aso_1161329 | 8 | 9 | 10.428518 | 14.204306 | 0.0e+00 | 0.0000001 | *** | 1215.6400 | SCRaso , MSH3aso_1161329 |
| MSH3_counts_normalised | MSH3aso_1161149 | MSH3aso_1161173 | 9 | 9 | -7.666037 | 13.398977 | 3.0e-06 | 0.0000042 | *** | 1297.9067 | MSH3aso_1161149, MSH3aso_1161173 |
| MSH3_counts_normalised | MSH3aso_1161149 | MSH3aso_1161329 | 9 | 9 | -3.886817 | 8.178619 | 4.0e-03 | 0.0050000 | ** | 1380.1733 | MSH3aso_1161149, MSH3aso_1161329 |
| MSH3_counts_normalised | MSH3aso_1161173 | MSH3aso_1161329 | 9 | 9 | -2.343382 | 8.459657 | 4.6e-02 | 0.0511111 | ns | 1462.4400 | MSH3aso_1161173, MSH3aso_1161329 |
# ANOVA and corrected posthoc pairwise comparisons
my_aov <- aov(MSH3_counts_normalised ~ treatment, data = plot_data)
my_emmeans <- emmeans(my_aov, pairwise ~ treatment) # emmeans 'adjust' argument is giving odd results, with some p values getting more significant
my_pwc <- as.data.frame(my_emmeans$contrasts) %>%
dplyr::mutate(group1 = str_split_fixed(contrast, " - ", 2)[, 1],
group2 = str_split_fixed(contrast, " - ", 2)[, 2],
p.adj = p.adjust(p.value, method = padj_method),
p.adj.signif = symnum(p.adj, corr = FALSE, na = FALSE, legend = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns"))) %>%
relocate(c(group1, group2), .after = "contrast") %>%
relocate(c(p.adj, p.adj.signif), .after = "p.value")
if (hide_ns) {
my_pwc_plot <- my_pwc %>%
dplyr::filter(p.adj < 0.05)
} else {
my_pwc_plot <- my_pwc
}
y_max <- max(plot_data$MSH3_counts_normalised, na.rm = TRUE)
if (nrow(my_pwc_plot) > 0) {
my_pwc_plot <- my_pwc_plot %>%
dplyr::mutate(ymax = y_max) %>%
dplyr::mutate(y.position = ymax + seq(0.1 * unique(ymax), 0.5 * unique(ymax), length.out = n()))
} else {
my_pwc_plot[1,] <- NA
my_pwc_plot$y.position <- NA
}
# Display ANOVA posthoc
kbl(my_pwc, caption = "ANOVA with corrected posthoc pairwise comparisons") %>%
kable_styling(bootstrap_options = "striped")| contrast | group1 | group2 | estimate | SE | df | t.ratio | p.value | p.adj | p.adj.signif |
|---|---|---|---|---|---|---|---|---|---|
| Vehicle - SCRaso | Vehicle | SCRaso | 68.97222 | 25.67849 | 39 | 2.685992 | 0.0744751 | 0.0930938 | ns |
| Vehicle - MSH3aso_1161149 | Vehicle | MSH3aso_1161149 | 504.77778 | 24.91179 | 39 | 20.262605 | 0.0000000 | 0.0000000 | *** |
| Vehicle - MSH3aso_1161173 | Vehicle | MSH3aso_1161173 | 471.55556 | 24.91179 | 39 | 18.929010 | 0.0000000 | 0.0000000 | *** |
| Vehicle - MSH3aso_1161329 | Vehicle | MSH3aso_1161329 | 420.00000 | 24.91179 | 39 | 16.859486 | 0.0000000 | 0.0000000 | *** |
| SCRaso - MSH3aso_1161149 | SCRaso | MSH3aso_1161149 | 435.80556 | 25.67849 | 39 | 16.971622 | 0.0000000 | 0.0000000 | *** |
| SCRaso - MSH3aso_1161173 | SCRaso | MSH3aso_1161173 | 402.58333 | 25.67849 | 39 | 15.677845 | 0.0000000 | 0.0000000 | *** |
| SCRaso - MSH3aso_1161329 | SCRaso | MSH3aso_1161329 | 351.02778 | 25.67849 | 39 | 13.670112 | 0.0000000 | 0.0000000 | *** |
| MSH3aso_1161149 - MSH3aso_1161173 | MSH3aso_1161149 | MSH3aso_1161173 | -33.22222 | 24.91179 | 39 | -1.333594 | 0.6723546 | 0.6723546 | ns |
| MSH3aso_1161149 - MSH3aso_1161329 | MSH3aso_1161149 | MSH3aso_1161329 | -84.77778 | 24.91179 | 39 | -3.403119 | 0.0127230 | 0.0181758 |
|
| MSH3aso_1161173 - MSH3aso_1161329 | MSH3aso_1161173 | MSH3aso_1161329 | -51.55556 | 24.91179 | 39 | -2.069524 | 0.2536191 | 0.2817990 | ns |
# Plot
my_plot <-
ggplot(plot_data, aes(x = treatment, y = MSH3_counts_normalised)) +
# stat_summary(fun = mean, geom = "bar", colour = "black", alpha = 0.8) +
stat_summary(fun = mean, geom = "bar", colour = "black") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
geom_jitter(width = 0.2, height = 0, size = point_size, alpha = 0.6) +
stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
# geom_hline(yintercept = 0, colour = "darkgrey") +
scale_x_discrete(labels = treatment_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Treatment",
y = "Normalised MSH3 counts") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, ".", output_format)),
plot = my_plot,
height = 6, width = 6)
}
# Plot in colour
my_plot <-
ggplot(plot_data, aes(x = treatment, y = MSH3_counts_normalised)) +
stat_summary(aes(fill = treatment),
# fun = mean, geom = "bar", colour = "black", alpha = 0.8) +
fun = mean, geom = "bar", colour = "black") +
stat_summary(fun.data = mean_se, geom = "errorbar", width = 0.2) +
geom_jitter(aes(fill = treatment),
shape = 21, width = 0.2, height = 0, size = point_size, alpha = 0.8) +
stat_pvalue_manual(my_pwc_plot, hide.ns = hide_ns, label = "p.adj.signif", tip.length = 0.01) +
# geom_hline(yintercept = 0, colour = "darkgrey") +
scale_fill_manual(values = treatment_colour,
labels = treatment_rename) +
scale_x_discrete(labels = treatment_rename,
guide = guide_axis(angle = 45)) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "Treatment",
y = "Normalised MSH3 counts") +
theme_bw() +
theme(text = element_text(size = element_text_size),
legend.position = "none")
my_plot
MSH3 protein vs instability
Format data
## [1] "MSH3_western_instab_comparison"# Make results directory
dir.create(file.path(out_dir, experiment_id),
recursive = TRUE)
# Format data
plot_data <- MSH3aso_data[[experiment_id]] %>%
dplyr::mutate(clone = "") %>%
relocate(clone, .after = "genotype")
# Add MSH3ko protein level
MSH3ko_protein <- data.frame("experiment_id" = "MSH3ko_CRISPRwt",
"differentiation_id" = c("NI1908", "NI1708", "NI1808"),
"genotype" = "MSH3ko",
"clone" = c("Cl27", "Cl37", "Cl26"),
"condition_eb" = NA,
"treatment" = "untreated",
"dose_uM" = NA,
"treatment_dose" = "untreated",
"group_juliet" = "MSH3ko_untreated",
"group_juliet_differentiations" = c("MSH3ko_CRISPRwt-NI1908-Cl27-MSH3ko",
"MSH3ko_CRISPRwt-NI1708-Cl37 -MSH3ko",
"MSH3ko_CRISPRwt-NI1808-Cl26-MSH3ko"),
"MSH3_level_normalised" = 0)
plot_data <- rbind(plot_data, MSH3ko_protein) %>%
dplyr::mutate(predvar2 = paste0(experiment_id, "-", differentiation_id, "-", clone, "-", genotype))
# Add instability rate
plot_data <- plot_data %>%
left_join(juliet_slopes %>%
dplyr::mutate(predvar2 = str_replace_all(predvar2, " ", "")),
by = join_by(predvar2, group_juliet == predvar1)) %>%
dplyr::rename(ii_slope = slope)
# Add annotation
plot_data <- plot_data %>%
left_join(juliet_annotation,
by = join_by(experiment_id, group_juliet)) %>%
dplyr::mutate(group_juliet = fct_relevel(group_juliet, unique(juliet_annotation$group_juliet)))Fit models
# Plot using each formula
my_fits <- lapply(model_names, plot_model_function,
models = my_models,
data = plot_data,
xvar = "MSH3_level_normalised",
xlabel = "MSH3 protein level (%)",
yvar = "ii_slope",
ylabel = "CAG expansion rate (Change in instability index / week)",
group_var = NA,
group_label = NA,
log_increment = log_increment,
plot_se = TRUE)
# Plot model fits
my_plots <- lapply(my_fits, function(my_fit) {
my_fit$my_plot
})
wrap_plots(my_plots) +
plot_annotation(title = "Regression models",
theme = theme(plot.title = element_text(hjust = 0.5, face = "bold")))
# r2 table
my_r2 <- lapply(my_fits, function(my_fit) {
my_fit[["r2"]]
})
r2_table <- data.frame(model_name = names(my_r2),
r2 = unlist(my_r2)) %>%
tibble::rownames_to_column("model_number") %>%
arrange(-r2)
# Display r2 table
kbl(r2_table, caption = "R-squared") %>%
kable_styling(bootstrap_options = "striped")| model_number | model_name | r2 |
|---|---|---|
| 4 | poly3 | 0.8377280 |
| 3 | poly2 | 0.8335948 |
| 1 | linear | 0.7820496 |
| 2 | log | 0.5404600 |
| 5 | exponential | 0.4965232 |
| 6 | exponential | 0.4965232 |
IC50
Dose response model
# Fit dose-response model
my_drm <- drm(ii_slope ~ MSH3_level_normalised, data = plot_data, fct = LL.4())
# Extract IC50 value
ic50 <- ED(my_drm, 50)##
## Estimated effective doses
##
## Estimate Std. Error
## e:1:50 32.290 16.092ic50_estimate <- ic50[[1]]
ic50_se <- ic50[[2]]
ic50_upper <- ic50_estimate + 1.96 * ic50_se
ic50_lower <- ic50_estimate - 1.96 * ic50_se
# Find max value for the response variable using the dose response model
max_response_var_drm <- predict(my_drm, newdata = data.frame(MSH3_level_normalised = 100))
# Plot drm model
plot(my_drm, main = "Dose-response curve using drc")
abline(h = 0.5 * max_response_var_drm, col = "blue", lty = 2)
abline(v = ic50_estimate, col = "blue", lty = 2)
# Table of IC50
ic50_table <- data.frame(method = "drc",
max_ii_slope = as.numeric(max_response_var_drm),
IC50 = ic50_estimate,
ci_lower = ic50_lower,
ci_upper = ic50_upper)
kbl(ic50_table, caption = "IC50") %>%
kable_styling(bootstrap_options = "striped")| method | max_ii_slope | IC50 | ci_lower | ci_upper |
|---|---|---|---|---|
| drc | 0.105816 | 32.29007 | 0.7495401 | 63.83061 |
Interpolation using a dose response model
# Use predicted data from drm curve to calculate the drm curve IC50 by interpolation
# predicted_values_drm <- data.frame(MSH3_level_normalised = seq(0, max(plot_data$MSH3_level_normalised), length.out = 1000))
predicted_values_drm <- data.frame(MSH3_level_normalised = seq(0, max(plot_data$MSH3_level_normalised), length.out = nrow(plot_data)))
predictions <- as.data.frame(predict(my_drm, newdata = predicted_values_drm, interval = "confidence")) %>%
dplyr::rename(ii_slope = Prediction,
ci_lower = Lower,
ci_upper = Upper)
predicted_values_drm <- cbind(predicted_values_drm, predictions)
# plot(ii_slope ~ MSH3_level_normalised, data = predicted_values_drm)
# Calculate IC50 using interpolation with the dose response model
ic50_interpolation_drm <- ed_interpolation(data = predicted_values_drm,
response_var = "ii_slope",
pred_var = "MSH3_level_normalised",
formula = as.formula("ii_slope ~ MSH3_level_normalised"),
model_function = "drm",
response_level = 0.5 * max_response_var_drm,
n_bootstrap = 1000,
confidence_level = 0.95)
# Calculate IC100 (response variable 0% of its maximum value)
ic100_interpolation_drm <- ed_interpolation(data = predicted_values_drm,
response_var = "ii_slope",
pred_var = "MSH3_level_normalised",
formula = as.formula("ii_slope ~ MSH3_level_normalised"),
model_function = "drm",
response_level = 0,
n_bootstrap = 1000,
confidence_level = 0.95)
# Plot linear model
plot(ii_slope ~ MSH3_level_normalised,
data = ic50_interpolation_drm$predicted_values,
main = "Dose-response curve using interpolation")
abline(h = 0.5 * max_response_var_drm, col = "blue", lty = 2)
abline(v = ic50_interpolation_drm$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation_drm$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation_drm$ci_upper, col = "lightblue", lty = 2)
abline(h = 0, col = "red", lty = 2)
abline(v = ic100_interpolation_drm$ed, col = "red", lty = 2)
abline(v = ic100_interpolation_drm$ci_lower, col = "pink", lty = 2)
abline(v = ic100_interpolation_drm$ci_upper, col = "pink", lty = 2)
# Table of IC50
ed_table <- data.frame(method = "drc_interpolation",
max_ii_slope = as.numeric(max_response_var_drm),
ed_level = c(50, 100),
ED = c(ic50_interpolation_drm$ed, ic100_interpolation_drm$ed),
ci_lower = c(ic50_interpolation_drm$ci_lower, ic100_interpolation_drm$ci_lower),
ci_upper = c(ic50_interpolation_drm$ci_upper, ic100_interpolation_drm$ci_upper))
kbl(ed_table, caption = "Effective doses (ED)") %>%
kable_styling(bootstrap_options = "striped")| method | max_ii_slope | ed_level | ED | ci_lower | ci_upper |
|---|---|---|---|---|---|
| drc_interpolation | 0.105816 | 50 | 36.81592 | 36.81562 | 36.82542 |
| drc_interpolation | 0.105816 | 100 | 18.13204 | 18.07901 | 18.13435 |
Interpolation using a linear model
# Calculate IC50 by interpolation
my_lm <- lm(ii_slope ~ MSH3_level_normalised, data = plot_data)
max_response_var_lm <- predict(my_lm, newdata = data.frame(MSH3_level_normalised = 100))
ic50_interpolation_lm <- ed_interpolation(data = plot_data,
response_var = "ii_slope",
pred_var = "MSH3_level_normalised",
formula = as.formula("ii_slope ~ MSH3_level_normalised"),
model_function = "lm",
response_level = 0.5 * max_response_var_lm,
n_bootstrap = 1000,
confidence_level = 0.95)
# Calculate IC100 (response variable 0% of its maximum value)
ic100_interpolation_lm <- ed_interpolation(data = plot_data,
response_var = "ii_slope",
pred_var = "MSH3_level_normalised",
formula = as.formula("ii_slope ~ MSH3_level_normalised"),
model_function = "lm",
response_level = 0,
n_bootstrap = 1000,
confidence_level = 0.95)
# Plot linear model
plot(ii_slope ~ MSH3_level_normalised,
data = ic50_interpolation_lm$predicted_values,
main = "Dose-response curve using interpolation")
abline(h = 0.5 * max_response_var_lm, col = "blue", lty = 2)
abline(v = ic50_interpolation_lm$ed, col = "blue", lty = 2)
abline(v = ic50_interpolation_lm$ci_lower, col = "lightblue", lty = 2)
abline(v = ic50_interpolation_lm$ci_upper, col = "lightblue", lty = 2)
abline(h = 0, col = "red", lty = 2)
abline(v = ic100_interpolation_lm$ed, col = "red", lty = 2)
abline(v = ic100_interpolation_lm$ci_lower, col = "pink", lty = 2)
abline(v = ic100_interpolation_lm$ci_upper, col = "pink", lty = 2)
# Table of effective doses
ed_table <- data.frame(method = "lm_interpolation",
max_ii_slope = as.numeric(max_response_var_lm),
ed_level = c(50, 100),
ED = c(ic50_interpolation_lm$ed, ic100_interpolation_lm$ed),
ci_lower = c(ic50_interpolation_lm$ci_lower, ic100_interpolation_lm$ci_lower),
ci_upper = c(ic50_interpolation_lm$ci_upper, ic100_interpolation_lm$ci_upper))
kbl(ed_table, caption = "Effective doses (ED)") %>%
kable_styling(bootstrap_options = "striped")| method | max_ii_slope | ed_level | ED | ci_lower | ci_upper |
|---|---|---|---|---|---|
| lm_interpolation | 0.1228149 | 50 | 58.71070 | 39.227615 | 72.26981 |
| lm_interpolation | 0.1228149 | 100 | 17.42139 | 4.648074 | 27.23633 |
Dose response curve
# Group colours
my_colours <- juliet_annotation %>%
dplyr::select(group_juliet, group_colour) %>%
distinct()
my_colours <- setNames(my_colours$group_colour, my_colours$group_juliet)
# Group rename
my_rename <- juliet_annotation %>%
dplyr::select(group_juliet, group_rename) %>%
distinct()
my_rename <- setNames(my_rename$group_rename, my_rename$group_juliet)
# Plot
my_plot <-
ggplot(plot_data, aes(x = MSH3_level_normalised, y = ii_slope)) +
geom_point(aes(fill = group_juliet),
size = point_size, shape = 21, colour = "black", alpha = 0.8) +
geom_line(data = ic50_interpolation_drm$predicted_values,
aes(x = MSH3_level_normalised, y = ii_slope),
colour = "black", size = 1) +
geom_segment(aes(x = -Inf, xend = ic50_interpolation_drm$ed,
y = 0.5 * max_response_var_drm, yend = 0.5 * max_response_var_drm),
colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_segment(aes(x = ic50_interpolation_drm$ed, xend = ic50_interpolation_drm$ed,
y = -Inf, yend = 0.5 * max_response_var_drm),
colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_rect(aes(xmin = ic50_interpolation_drm$ci_lower, xmax = ic50_interpolation_drm$ci_upper,
ymin = -Inf, ymax = 0.5 * max_response_var_drm),
fill = "blue", alpha = 0.01, colour = NA) +
geom_segment(aes(x = -Inf, xend = ic100_interpolation_drm$ed,
y = 0, yend = 0),
colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_segment(aes(x = ic100_interpolation_drm$ed, xend = ic100_interpolation_drm$ed,
y = -Inf, yend = 0),
colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_rect(aes(xmin = ic100_interpolation_drm$ci_lower, xmax = ic100_interpolation_drm$ci_upper,
ymin = -Inf, ymax = 0),
fill = "red", alpha = 0.01, colour = NA) +
scale_fill_manual(values = my_colours,
labels = my_rename) +
scale_x_continuous(breaks = pretty_breaks()) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 protein level (%)",
y = "CAG expansion rate\n(change in instability index / week)",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot
# Export
for (output_format in output_formats) {
ggsave(filename = file.path(out_dir, experiment_id, paste0(experiment_id, "_drm.", output_format)),
plot = my_plot,
height = 6, width = 9)
}
# Add confidence interval
my_plot <-
my_plot +
geom_ribbon(data = predicted_values_drm,
aes(ymin = ci_lower, ymax = ci_upper),
fill = "darkgrey", alpha = 0.3)
my_plot
Linear model
my_plot <-
ggplot(plot_data, aes(x = MSH3_level_normalised, y = ii_slope)) +
geom_point(aes(fill = group_juliet),
size = point_size, shape = 21, colour = "black", alpha = 0.8) +
stat_smooth(method = lm,
formula = y ~ x,
fullrange = TRUE, se = TRUE, colour = "black", fill = "darkgrey", alpha = 0.3) +
geom_segment(aes(x = -Inf, xend = ic50_interpolation_lm$ed,
y = 0.5 * max_response_var_lm, yend = 0.5 * max_response_var_lm),
colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_segment(aes(x = ic50_interpolation_lm$ed, xend = ic50_interpolation_lm$ed,
y = -Inf, yend = 0.5 * max_response_var_lm),
colour = "blue", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_rect(aes(xmin = ic50_interpolation_lm$ci_lower, xmax = ic50_interpolation_lm$ci_upper,
ymin = -Inf, ymax = 0.5 * max_response_var_lm),
fill = "blue", alpha = 0.01, colour = NA) +
geom_segment(aes(x = -Inf, xend = ic100_interpolation_lm$ed,
y = 0, yend = 0),
colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_segment(aes(x = ic100_interpolation_lm$ed, xend = ic100_interpolation_lm$ed,
y = -Inf, yend = 0),
colour = "red", linetype = "dashed", size = 0.5, alpha = 0.5) +
geom_rect(aes(xmin = ic100_interpolation_lm$ci_lower, xmax = ic100_interpolation_lm$ci_upper,
ymin = -Inf, ymax = 0),
fill = "red", alpha = 0.01, colour = NA) +
scale_fill_manual(values = my_colours,
labels = my_rename) +
scale_x_continuous(breaks = pretty_breaks()) +
scale_y_continuous(breaks = pretty_breaks()) +
labs(x = "MSH3 protein level (%)",
y = "CAG expansion rate\n(change in instability index / week)",
fill = "Treatment") +
theme_bw() +
theme(text = element_text(size = element_text_size))
my_plot